############################################################
#
#   Sovereign Debt Ratchet and Welfare Destruction
#   July 2022
#
#############################################################

#================
# Basic Functions
#================

tic <- function(gcFirst = TRUE, type=c("elapsed", "user.self", "sys.self"))
{
  type <- match.arg(type)
  assign(".type", type, envir=baseenv())
  if(gcFirst) gc(FALSE)
  tic <- proc.time()[type]         
  assign(".tic", tic, envir=baseenv())
  invisible(tic)
}

toc <- function()
{
  type <- get(".type", envir=baseenv())
  toc <- proc.time()[type]
  tic <- get(".tic", envir=baseenv())
  print(toc - tic)
  invisible(toc)
}

FB_fn <- function(x,q,M,l,u){
  # Fischer-Burmeister function
  # output of the FB function: psi,phi,J
  # psi is a scalar, phi is an n x 1 array and J is an n x n matrix
  # inputs are as follows:
  # x: n x 1 vector
  # M: n x n matrix
  # q: n x 1 vector
  # l: n x 1 vector
  # u: n x 1 vector
  
  n     <- length(x)
  Zl    <- (l > -Inf)  & (u == Inf)  # this is the usual case when l and u are not specified
  Zu    <- (l == -Inf) & (u < Inf)   # check for non-trivial upper boundaries and no lower boundaries
  Zlu   <- (l > -Inf)  & (u < Inf)   # check for non-trivial upper and lower boundaries
  Zf    <- (l== -Inf)  & (u == Inf)
  
  a     <- x
  b     <- c(array(M %*% x + q))
  
  a[Zl] <- x[Zl]-l[Zl]
  a[Zu] <- u[Zu]-x[Zu]
  b[Zu] <- -b[Zu]
  
  if (sum(Zlu) >= 1){
    nt     <- sum(Zlu)
    at     <- u[Zlu]-x[Zlu]
    bt     <- -b[Zlu]
    st     <- sqrt(at^2 + bt^2)
    a[Zlu] <- x[Zlu]-l[Zlu]
    b[Zlu] <- st - at - bt
  }
  
  s        <- c(sqrt(a^2 + b^2))
  phi      <- s - a - b
  phi[Zu]  <- -phi[Zu]
  phi[Zf]  <- -b[Zf]
  psi      <- 0.5*sum(phi^2)
  
  if (sum(Zlu) >= 1){
    M1       <- Matrix(0, nrow = nt, ncol = n, sparse = TRUE)
    M2       <- Matrix(0, nrow = nt, ncol = nt, sparse = TRUE)
    loc1     <- cbind(1:nt, which(Zlu))
    M1[loc1] <- at/st - rep(1,nt)
    loc2     <- cbind(seq(nt), seq(nt))
    M2[loc2] <- bt/st - rep(1,nt)
    M[Zlu,]  <- -M1 - M2 %*% M[Zlu,]
  }
  da       <- a/s-rep(1,n)
  db       <- b/s-rep(1,n)
  da[Zf]   <- 0
  db[Zf]   <- -1
  J1       <- Matrix(0, nrow = n, ncol = n, sparse = TRUE)
  J1[cbind(seq(n),seq(n))] <- da
  J2       <- Matrix(0, nrow = n, ncol = n, sparse = TRUE)
  J2[cbind(seq(n),seq(n))] <- db
  J        <- J1 + J2 %*% M
  
  output         <- list(psi,phi,J)
  names(output)  <- c('psi','phi','J')
  return(output)
}

LCP_fn <- function(M,q,l,u,x0,display){
  # LCP Solve the Linear Complementarity Problem.
  # 
  #  USAGE
  #    x = LCP(M,q) solves the LCP
  # 
  #            x >= 0
  #       Mx + q >= 0 
  #    x'(Mx + q) = 0  
  # 
  #    x = LCP(M,q,l,u) solves the generalized LCP (a.k.a MCP)
  # 
  #    l < x < u   =>   Mx + q = 0
  #        x = u   =>   Mx + q < 0
  #    l = x       =>   Mx + q > 0
  # 
  #    x = LCP(M,q,l,u,x0,display) allows the optional initial value 'x0' and
  #    a binary flag 'display' which controls the display of iteration data.
  # 
  #    Parameters:
  #    tol       -   Termination criterion. Exit when 0.5*phi(x)'*phi(x) < tol.
  #    mu        -   Initial value of Levenberg-Marquardt mu coefficient.
  #    mu_step   -   Coefficient by which mu is multiplied / divided.
  #    mu_min    -   Value below which mu is set to zero (pure Gauss-Newton).
  #    max_iter  -   Maximum number of (succesful) Levenberg-Marquardt steps.
  #    b_tol     -   Tolerance of degenerate complementarity: Dimensions where
  #                  max( min(abs(x-l),abs(u-x)) , abs(phi(x)) ) < b_tol
  #                  are clamped to the nearest constraint and removed from
  #                  the linear system.
  #    
  #  ALGORITHM
  #    This function implements the semismooth algorithm as described in [1],
  #    with a least-squares minimization of the Fischer-Burmeister function using
  #    a Levenberg-Marquardt trust-region scheme with mu-control as in [2].
  # 
  #    [1] A. Fischer, A Newton-Type Method for Positive-Semidefinite Linear
  #    Complementarity Problems, Journal of Optimization Theory and
  #    Applications: Vol. 86, No. 3, pp. 585-608, 1995.
  # 
  #    [2] M. S. Bazarraa, H. D. Sherali, and C. M. Shetty, Nonlinear
  #    Programming: Theory and Algorithms. John Wiley and Sons, 1993.
  # 
  #    Copyright (c) 2008, Yuval Tassa
  #    tassa at alice dot huji dot ac dot il
  
  tol            <- 1.0e-12
  mu             <- 1e-3
  mu_step        <- 5
  mu_min         <- 1e-5
  max_iter       <- 20
  b_tol          <- 1e-6
  n              <- dim(M)[1]
  
  if ((nargs()<3)|(length(l)==0)){
    l <- rep(0,n)
  } else if ((nargs()<4)|(length(u)==0)){
    u <- rep(Inf,n)
  } else if ((nargs()<5)|(length(x0)==0)){
    x0 <- pmin(pmax(rep(1,n),l),u) 
  } else if (nargs()<6){
    display <- FALSE
  }
  
  # create matrix lu
  lu <- cbind(l,u)
  x  <- x0
  
  # compute psi,phi,J
  FB_out1 <- FB_fn(x,q,M,l,u)
  psi     <- FB_out1$psi       # this is a scalar
  phi     <- FB_out1$phi       # this is an n x 1 vector
  J       <- FB_out1$J         # this is an n x n matrix
  
  new_x  <- TRUE
  iter   <- 1
  while ((iter < max_iter)&(psi > tol)){
    
    if (new_x){
      # find the minimum of abs(x-l) and abs(u-x), row by row
      mlu            <- pmin(abs(x-l),abs(u-x))
      # find the index of the minimum of abs(x-l) and abs(u-x), row by row
      ilu            <- (mlu==abs(x-l))*1 + (1 - (mlu==abs(x-l)))*2
      bad            <- pmax(abs(phi),mlu) < b_tol
      psi            <- psi - 0.5*sum(phi[bad]^2)
      J              <- J[!bad,!bad]
      phi            <- phi[!bad] 
      new_x          <- FALSE
      nx             <- x
      nx[bad]        <- lu[which(bad)+(ilu[bad]-1)*n]
    }
    
    H              <- t(J) %*% J + mu*Diagonal(sum(!bad))
    Jphi           <- c(array(t(J) %*% phi))
    d              <- - c(array(solve(H,Jphi)))
    nx[!bad]       <- x[!bad] + d
    
    # compute npsi,nphi,nJ
    FB_out2 <- FB_fn(nx,q,M,l,u)
    npsi    <- FB_out2$psi
    nphi    <- FB_out2$phi
    nJ      <- FB_out2$J
    
    # calculate actual reduction / expected reduction
    r       <- - (psi - npsi)  / (sum(Jphi * d) + 0.5*c(array(t(d) %*% H %*% d)))  
    
    # if small reduction, increase mu
    if (r < 0.3) mu  <- max(mu*mu_step,mu_min)
    
    if (r > 0) {
      #some reduction, accept nx
      x     <- nx
      psi   <- npsi
      phi   <- nphi
      J     <- nJ
      new_x <- TRUE
      # if large reduction, decrease mu
      if (r > 0.8) mu <- mu/mu_step * (mu > mu_min)
    }
    
    if (display) print(paste('iter = ',iter,' -- psi = ',psi,' -- r = ',r,' -- mu = ',mu,sep=''))
    
    iter <- iter + 1
  }
  
  # return output
  x <- pmin(pmax(x,l),u)
  return(x)
}

#==================================
# Functions for No-Commitment Model
#==================================

no_commit_root_fn <- function(param){
  # compute root xi for utility function without commitment
  delta   <- param["delta"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  
  # compute root of characteristic polynomial
  a       <- 0.5*sigma^2
  b       <- - (m + mu + 0.5*sigma^2)
  c       <- - (delta - mu)
  discrim <- b^2 - 4*a*c
  xi      <- (- b + sqrt(discrim))/(2*a)
  
  return(xi)
}

no_commit_root2_fn <- function(param){
  # compute roots xi_1 and xi_2 for utility function without commitment
  # return the 2 roots
  delta   <- param["delta"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  
  # compute root of characteristic polynomial
  a       <- 0.5*sigma^2
  b       <- - (m + mu + 0.5*sigma^2)
  c       <- - (delta - mu)
  discrim <- b^2 - 4*a*c
  xi_1    <- (- b - sqrt(discrim))/(2*a)
  xi_2    <- (- b + sqrt(discrim))/(2*a)
  xi      <- c(xi_1,xi_2)
  names(xi) <- c("xi_1", "xi_2")
  return(xi)
}

no_commit_attraction_pt_fn <- function(param){
  # compute attraction point for case without commitment
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  nu      <- param["nu"]
  
  # compute root of characteristic polynomial and default boundary
  xi      <- no_commit_root_fn(param)
  x_bar   <- no_commit_cutoff_fn(param)
  
  # compute point at which drift = 0
  att_pt  <- x_bar * ((1-alpha*theta)/(1-alpha*theta^xi)*((xi-1)*(m+nu*sigma)/
                                   (delta-r)+1))^(1/(1 - xi))
  
  return(att_pt)
}

no_commit_zero_drift_pt_fn <- function(param){
  # compute zero drift point mu(x)=0 for case without commitment
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  nu      <- param["nu"]
  
  # compute root of characteristic polynomial and default boundary
  xi      <- no_commit_root_fn(param)
  x_bar   <- no_commit_cutoff_fn(param)
  
  # compute point at which drift = 0
  x0      <- x_bar * ((1-alpha*theta)/(1-alpha*theta^xi)*((xi-1)*(m+nu*sigma+mu-sigma^2)/
                                                            (delta-r)+1))^(1/(1 - xi))
  
  return(x0)
}

no_commit_zero_ca_pt_fn <- function(param){
  # compute point c(x)=1 for case without commitment
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  nu      <- param["nu"]
  
  # compute default point
  x_bar   <- no_commit_cutoff_fn(param)
  
  # create implicit function
  imp_fn <- function(x_guess){
    return(no_commit_c_fn(param, x_guess)-1)
  }
  
  # solve implicit function
  x_zero_ca <- uniroot(f = imp_fn, interval = c(.1,x_bar))$root
  
  return(x_zero_ca)
}

no_commit_attraction_pt_fn2 <- function(param){
  # compute attraction point for case without commitment
  # state variable is now y
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  nu      <- param["nu"]
  
  # compute root of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  x_bar   <- no_commit_cutoff_fn(param)
  y_bar   <- 1 / x_bar
  
  # compute point at which drift = 0
  att_pt  <- y_bar * ((1-alpha*theta)/(1-alpha*theta^xi)*((xi-1)*(m+mu+nu*sigma)/
                                                            (delta-r)+1))^(1/(xi-1))
  
  return(att_pt)
}

no_commit_cutoff_fn <- function(param){
  # compute optimal default cutoff for case without commitment
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  
  # compute root of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  x_bar   <- xi/(xi-1) * ((delta + m)/(kappa + m)) * ((1 - alpha)/(1 - alpha * theta)) / (delta - mu)
  
  return(x_bar)
}

no_commit_debt_px_fn <- function(param, x){
  # compute debt price for case without commitment
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  eta     <- param["eta"]
  
  # compute root of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  
  # compute default boundary
  x_bar   <- no_commit_cutoff_fn(param)
  
  # compute d_0, risk free debt price
  d0      <- (kappa + m)/(delta + m)
  
  # compute debt price
  debt    <- (d0/(1-eta)) * (1 - (1 - alpha*theta)/(1 - alpha*(theta^xi)) * (x/x_bar)^(xi - 1))
  
  # take care of situations where x > x_bar -- compute number if consecutive restructurings
  n       <- pmax(0,1 + floor((log(x_bar) - log(x))/log(theta)))
  
  # compute debt price in case x > x_bar
  debt_def<- ((alpha*theta)^n) * (d0/(1-eta)) * (1 - (1 - alpha*theta)/(1 - alpha*(theta^xi)) * ((theta^n)*x/x_bar)^(xi - 1))
  
  # return debt price
  out <- debt*(n<=0) + debt_def*(n>0)
  return(out)
}

no_commit_value_fn <- function(param, x){
  # compute value function for case without commitment
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  
  # compute root of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  # compute default boundary
  x_bar   <- no_commit_cutoff_fn(param)
  # compute value function
  value   <- (1 / (delta - mu) * 
                    (1 - (1 - alpha)/(1 - alpha*(theta^xi)) * (x/x_bar)^(xi)) - 
                    x * (kappa + m)/(delta + m) * 
                    (1 - (1 - alpha*theta)/(1 - alpha*(theta^xi)) * (x/x_bar)^(xi - 1)))
  # take care of situations where x > x_bar -- compute number if consecutive restructurings
  n       <- pmax(0,1 + floor((log(x_bar) - log(x))/log(theta)))
  # compute value in case x > x_bar
  val_def <- (alpha^n) * (1 / (delta - mu) * 
                            (1 - (1 - alpha)/(1 - alpha*(theta^xi)) * ((theta^n)*x/x_bar)^(xi)) - 
                            (theta^n)*x * (kappa + m)/(delta + m) * 
                            (1 - (1 - alpha*theta)/(1 - alpha*(theta^xi)) * ((theta^n)*x/x_bar)^(xi - 1)))
  # return value
  val <- value*(n<=0) + val_def*(n>0)
  return(val)
}

no_commit_g_fn <- function(param, x){
  # compute issuance policy for case without commitment
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  nu      <- param["nu"]
  eta     <- param["eta"]
  
  # compute root of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  
  # compute default boundary
  x_bar   <- no_commit_cutoff_fn(param)
  
  # compute issuance policy
  g       <- x * (delta - r) / (xi - 1) * (((1 - alpha*(theta^xi))/(1 - alpha*theta)) * 
       (1 - eta*(delta+m)/(1-eta)/(delta-r)) * (x_bar/x)^(xi-1) - 1) - x * nu * sigma
  
  return(g)
}

no_commit_c_fn <- function(param, x){
  # compute consumption-to-income policy for case without commitment
  delta      <- param["delta"]
  kappa      <- param["kappa"]
  m          <- param["m"]
  mu         <- param["mu"]
  sigma      <- param["sigma"]
  alpha      <- param["alpha"]
  theta      <- param["theta"]
  r          <- param["r"]
  nu         <- param["nu"]
  eta        <- param["eta"]
  deadweight <- param["deadweight_flag"]
  
  # dead-weight loss in case not specified
  if (is.na(deadweight)) deadweight <- 0
  
  # compute root of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  
  # compute risk free debt price
  d0      <- (kappa+m)/(delta+m)
  
  # compute default boundary
  x_bar   <- no_commit_cutoff_fn(param)
  
  # compute issuance policy  
  g       <- x * (delta - r) / (xi - 1) * (((1 - alpha*(theta^xi))/(1 - alpha*theta)) * 
       (1 - eta*(delta+m)/(1-eta)/(delta-r)) * (x_bar/x)^(xi-1) - 1) - x * nu * sigma
  
  # compute debt price
  debt    <- (d0/(1-eta)) * (1 - (1 - alpha*theta)/(1 - alpha*(theta^xi)) * (x/x_bar)^(xi - 1))
  
  # compute consumption policy
  c       <- 1 - (kappa+m)*x + debt*g - deadweight*eta*debt*g
  return(c)
}

no_commit_drift_fn <- function(param, x){
  # compute drift function for case without commitment
  # in this case we are using the state variable x := F / Y
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  nu      <- param["nu"]
  eta     <- param["eta"]
  
  # compute root of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  
  # compute default boundary
  x_bar   <- no_commit_cutoff_fn(param)
  
  # compute drift
  drift    <- (delta - r) * x / (xi - 1) * 
    ((1 - alpha*theta^xi)/(1 - alpha*theta) * (1 - eta*(delta+m)/(1-eta)/(delta-r)) * (x_bar/x)^(xi-1) - 1) - 
    x * (nu * sigma + m + mu - sigma^2)
  
  return(drift)
}

no_commit_drift_fn2 <- function(param, y){
  # compute drift function for case without commitment
  # in this case we are using the state variable y := Y / F
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  nu      <- param["nu"]
  eta     <- param["eta"]
  
  # compute root of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  
  # compute default boundary
  x_bar   <- no_commit_cutoff_fn(param)
  y_bar   <- 1 / x_bar
  
  # compute drift
  drift    <- y * (nu * sigma + m + mu) - 
    (delta - r) * (y / (xi - 1)) * ((1 - alpha*theta^xi)/(1 - alpha*theta) * (1 - eta*(delta+m)/(1-eta)/(delta-r)) * (y/y_bar)^(xi-1) - 1)
  
  return(drift)
}

no_commit_credit_spread_fn <- function(param, x){
  # compute credit spread function for case without commitment
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  nu      <- param["nu"]
  eta     <- param["eta"]
  
  # compute root of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  
  # compute default boundary
  x_bar   <- no_commit_cutoff_fn(param)
  
  # compute credit spread of exponentially amortizing bond
  spread    <- (m + delta)*(1 - eta) / 
    (1 - (1 - alpha*theta)/(1 - alpha*theta^xi) * 
       (x/x_bar)^(xi-1)) - (r + m)
  
  return(spread)
}

no_commit_risk_premium_fn <- function(param, x){
  # compute risk premium function for case without commitment
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  nu      <- param["nu"]
  eta     <- param["eta"]
  
  # compute root of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  
  # compute default boundary
  x_bar   <- no_commit_cutoff_fn(param)
  
  # compute risk-premium of exponentially amortizing bond
  risk_premium    <- sigma*nu*(xi-1)/
    (((1-alpha*theta^xi)/(1-alpha*theta))*(x_bar/x)^(xi-1)-1)
  
  return(risk_premium)
}

no_commit_cons_vol_fn <- function(param, x){
  # compute consumption vol for case without commitment
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  nu      <- param["nu"]
  eta     <- param["eta"]
  
  # compute root of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  
  # compute default boundary
  x_bar   <- no_commit_cutoff_fn(param)

  # compute consumption to income ratio
  c         <- no_commit_c_fn(param, x)
  
  # compute debt price and its derivative
  d0         <- (kappa+m)/(delta+m)
  debt       <- (d0/(1-eta)) * (1 - (1 - alpha*theta)/(1 - alpha*(theta^xi)) * (x/x_bar)^(xi - 1))
  debt_p     <- -(d0*(xi-1)/(1-eta)/x_bar) * (1 - alpha*theta)/(1 - alpha*(theta^xi)) * (x/x_bar)^(xi - 2)
  debt_pp    <- -(d0*(xi-1)*(xi-2)/(1-eta)/(x_bar^2)) * (1 - alpha*theta)/(1 - alpha*(theta^xi)) * (x/x_bar)^(xi - 3)
  
  # compute derivative of consumption to income ratio
  c_prime    <- - (kappa + m) - 2*(delta-r)*debt + (delta-r)*debt_pp*(debt/debt_p)^2 - nu*sigma*(debt + x*debt_p)
  
  # compute vol ratio
  cons_vol   <- sigma *(1 - x*c_prime/c)

  return(cons_vol)
}

no_commit_ergodic_dist_fn <- function(param){
  # compute ergodic distribution for case without commitment

  # load numerical parameters
  n           <- 2000
  dt          <- .25
  tol         <- .000001
  max_iter    <- 500
  dump_folder <- "C:/Users/ft.fi/Dropbox/Research/Sovereign Defaults and Commitment/Code/Dump"
  print_flag  <- F
  max_g       <- 1000  # to limit very intensities that make the matrix A no longer an intensity matrix
  max_cons    <- 1000  # to limit numerical instabilities arising from moment computation when we multiply very large numbers (consumption) by very small numbers (density)
  
  # load parameters
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  nu      <- param["nu"]
  eta     <- param["eta"]
  
  # compute equilibrium objects
  xi        <- no_commit_root_fn(param)
  x_bar     <- no_commit_cutoff_fn(param)
  
  # create x
  dx    <- x_bar / n
  dx2   <- dx^2
  
  # compute x_min
  # x_min <- (x_bar^(xi-1)*(delta-r)*(1 - alpha*theta^xi)*dx/sigma^2 / (xi-1) / (1 - alpha*theta))^(1/xi)
  x_min <- .01
  
  # create grid
  x     <- seq(from = x_min, to = x_bar, length.out = n)
  dx    <- x[2]-x[1]
  i_bar <- which.min(abs(x - x_bar))
  x_bar <- x[i_bar]
  
  # compute drift rate mu(x) and vol s(x) = sigma(x)^2
  g     <- pmin(no_commit_g_fn(param, x), max_g)
  drift <- g - (mu + m - sigma^2)*x
  
  # construct A matrix
  A      <- Matrix(0, nrow = n, ncol = n, sparse = TRUE)
  
  # diagonal elements
  locs   <- cbind(seq.int(n),seq.int(n))
  A[locs]<- -(abs(drift)/dx + (sigma*x)^2/dx2)
  
  # lower diagonal elements
  locs   <- cbind((2:n),(1:(n-1)))
  A[locs]<- -pmin(0,drift[-1])/dx + .5*(sigma*x[-1])^2/dx2
  
  # upper diagonal elements
  locs   <- cbind((1:(n-1)),(2:n))
  A[locs]<- pmax(0,drift[-n])/dx + .5*(sigma*x[-n])^2/dx2
  
  # modify intensity matrix with treatment of lower and higher boundary
  A[1,1] <- -sum(A[1,-1])
  A[n,n] <- -sum(A[n,-n])
  
  # check that rowsum = 0 for A_temp
  if (sum(abs(rowSums(A)))>.000001) print("Problem with A construction")

  # compute number of consecutive haircuts given state x. 0 if x < x_high for example
  haircut           <- rep(0,n)
  haircut[x>=x_bar] <- 1 + floor((log(x_bar)-log(x[x>=x_bar]))/log(theta))
  
  # compute states post default and restructuring (with possible multiple consecutive restructurings)
  i_next           <- seq(n)
  i_next[x>=x_bar] <- apply(as.matrix(i_bar:n), 1, FUN = function(u) return(which.min(abs(x - x[u]*theta^haircut[u]))))
  
  # compute intervention matrix M
  M <- Diagonal(n)
  
  # update intervention matrix for default -- reinject defaulted firms at lower debt-to-EBITDA ratio
  locs    <- cbind((i_bar:n),(i_bar:n))
  M[locs] <- 0
  locs    <- cbind((i_bar:n),i_next[i_bar:n])
  M[locs] <- 1
  
  # store transpose of intervention matrix and FK matrix
  tA      <- t(A)
  tM      <- t(M)
  
  # create/store B matrix
  B       <- (1/dt)*Diagonal(n) - tM %*% tA
  
  # initial guess stationary distribution -- assume everything centered at bar x / 2
  f_n                 <- rep(0,n)
  f_n[round(i_bar/2)] <- 1
  
  # initiate loop
  iter  <- 1
  err   <- 1
  err_t <- c(1)
  while ((err > tol) & (iter < max_iter)){
    # load previous period density
    f <- f_n
    
    if (print_flag) {
      # plot for verification purposes just in case
      setwd(dump_folder)
      png(paste("density_iter_",iter,".png",sep=""))
      plot(x = x, y = f_n, bty = 'n', type = 'l', xlim = c(0,x_bar), col = 'blue', xlab = 'debt-to-income', ylab = 'density')
      abline(v = x_bar, col = 'red', lty = 2)
      dev.off()
    }
    
    # update density from one loop to the next, and take into account interventions
    f   <- c(array(tM %*% f))
    
    # computing next step density using implicit scheme
    f_n <- c(array(solve(B,f/dt)))
    
    # normalize just in case
    f_n <- pmax(f_n,0)
    f_n <- f_n / sum(f_n)
    
    # normalize for distribution
    stat_dist <- f_n/dx
    
    # compute error
    err       <- max(abs(f_n-f)/dt)
    err_t     <- c(err_t,err)
    
    # update counter
    iter      <- iter + 1
  }
  
  # compute ergodic default rate
  lambda_d  <- .5*((sigma*x[i_bar-1])^2*stat_dist[i_bar-1]-(sigma*x[i_bar])^2*stat_dist[i_bar])/dx
  
  # compute stationary density
  stat_dist_fn  <- splinefun(x, stat_dist, method = "natural")
  
  # other moments we want to compute
  mom_stat             <- list()
  mom_stat[['x']]      <- x
  mom_stat[['z']]      <- (kappa+m)*x
  mom_stat[['spread']] <- no_commit_credit_spread_fn(param, x)
  mom_stat[['cons']]   <- pmin(max_cons,no_commit_c_fn(param, x))
  mom_stat[['ca']]     <- 1 - mom_stat[['cons']]
  mom_stat[['cons_vol']] <- no_commit_cons_vol_fn(param, x)
  mom_stat[['g']]      <- no_commit_g_fn(param, x)
  mom_stat[['v']]      <- no_commit_value_fn(param, x)
  mom_stat[['d']]      <- no_commit_debt_px_fn(param, x)

  # compute all statistics required
  moments           <- data.frame(matrix(0, nrow = length(mom_stat), ncol = 4))
  rownames(moments) <- names(mom_stat)
  colnames(moments) <- c('mean', 'stdev', 'skewness', 'kurtosis')
  for (mom in names(mom_stat)){
    moments[mom,'mean']     <- sum(f_n*mom_stat[[mom]])
    moments[mom,'stdev']    <- sqrt(sum(f_n*mom_stat[[mom]]^2) - moments[mom,'mean']^2)
    moments[mom,'skewness'] <- sum(f_n*((mom_stat[[mom]]-moments[mom,'mean'])/moments[mom,'stdev'])^3)
    moments[mom,'kurtosis'] <- sum(f_n*((mom_stat[[mom]]-moments[mom,'mean'])/moments[mom,'stdev'])^4)
  }
  
  # export all moments into csv file
  write.table(moments, file = paste("moments.csv",sep=""), eol = "\n", na = "NA", fileEncoding = "")
  
  # construct output
  output        <- list(stat_dist_fn, lambda_d, moments)
  names(output) <- c("density", "default_rate", "moments")
  return(output)
}

plot_fn <- function(param, plotting_param){
  # function that plots all useful plots
  
  # plotting parameters
  n              <- plotting_param['n']
  thickness      <- plotting_param['thick']
  font_size      <- plotting_param['font_size']
  axis_size      <- plotting_param['axis_size']
  
  # other parameters
  theta <- param['theta']
  alpha <- param['alpha']
  sigma <- param['sigma']
  
  # compute cutoff and other key cutoffs
  x_cutoff   <- no_commit_cutoff_fn(param)
  x_a        <- no_commit_attraction_pt_fn(param)
  x_0        <- no_commit_zero_drift_pt_fn(param)
  x_curr_act <- no_commit_zero_ca_pt_fn(param)
  x_reinject <- theta*x_cutoff
  
  # create grid
  if (2*param['delta']+param['m']+param['mu']-param['sigma']^2 > 0){
    # case where g(0)=+infty
    x_min    <- .5*theta*x_cutoff
  } else {
    # case where g(0)=0
    x_min    <- .0001
  }
  x_max    <- x_cutoff
  x        <- seq(from = x_min, to = x_max, length.out = n)
  
  # construct grid of output vectors
  out_list                <- list()
  out_list[['debt']]      <- no_commit_debt_px_fn(param, x)
  out_list[['value']]     <- no_commit_value_fn(param, x)
  out_list[['g']]         <- no_commit_g_fn(param, x)
  out_list[['drift']]     <- no_commit_drift_fn(param, x)
  out_list[['cons']]      <- no_commit_c_fn(param, x)
  out_list[['ca']]        <- 1 - no_commit_c_fn(param, x)
  out_list[['spread']]    <- no_commit_credit_spread_fn(param, x)
  out_list[['cons_vol']]  <- no_commit_cons_vol_fn(param, x)
  
  # stationary distribution
  dens_fn             <- no_commit_ergodic_dist_fn(param)[['density']]
  dens                <- dens_fn(x)
  
  # compute stationary mean
  x_mean              <- sum(dens*x)/sum(dens)
  
  # compute alpha * c(theta x_bar) = c post-default
  c_post_default      <- alpha*no_commit_c_fn(param, theta*x_cutoff)
  
  # reservation value at default
  res_val             <- alpha*no_commit_value_fn(param, theta*x)
  
  # debt price if commitment to never issue
  param_dummy          <- param
  param_dummy['gamma'] <- 0
  debt_commit          <- commit_debt_px_fn(param_dummy, x)
  spread_commit        <- commit_credit_spread_fn(param_dummy, x)
  risk_premium         <- no_commit_risk_premium_fn(param, x)
  
  # name plots
  legend <- c('debt price', 'value function', 'issuance policy', 'drift', 'consumption-to-income', 
              'current account-to-income', 'credit spread and debt risk-premium', 'consumption volatility')
  names(legend) <- names(out_list)

  # plot distribution on a standalone basis with key points
  pdf("stationary_density.pdf")
  par(mar=c(4, 4, 4, 4))
  plot(x = x, y = dens, ylim = c(0,1.4*max(dens)), type='l', xlab="debt-to-income x", 
       ylab="ergodic density", bty = 'n', 
       col="blue", main="", lwd = thickness, cex.lab = font_size, cex.axis = font_size)
  segments(x0 = x_mean, y0 = 0, x1 = x_mean, y1 = dens_fn(x_mean), lty = 2, col = 'purple', lwd = thickness)
  segments(x0 = x_0, y0 = 0, x1 = x_0, y1 = dens_fn(x_0), lty = 3, col = 'green', lwd = thickness)
  segments(x0 = x_a, y0 = 0, x1 = x_a, y1 = dens_fn(x_a), lty = 4, col = 'orange', lwd = thickness)
  segments(x0 = x_curr_act, y0 = 0, x1 = x_curr_act, y1 = dens_fn(x_curr_act), lty = 5, col = 'brown', lwd = thickness)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  legend("topright", col = c('purple','green','orange','brown'), bty="n", lty = 2:5, cex = font_size,
         legend = c('ergodic mean','drift(x) = 0','attraction point','current account(x) = 0'), lwd = thickness)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()

  # plot all functions of interest
  for (plot in names(out_list)){
    y_min <- min(out_list[[plot]], na.rm = TRUE)
    y_max <- max(out_list[[plot]], na.rm = TRUE)
    if ((plot == 'ca')|(plot == 'g')|(plot == 'drift')) {
      y_max <- min(1,y_max)
      y_min <- max(-1,y_min)
    }
    if (plot == 'cons') y_max = min(2,y_max)
    if (plot == 'cons_vol') y_min = 0
    if (plot == 'spread') {
      y_min = 0
      y_max = min(.05,y_max)
    }
    if (plot == 'debt') {
      y_min = 0
      y_max = 1
    }
    
    # setEPS()
    # postscript(paste(plot,".eps",sep=''))
    pdf(paste(plot,".pdf",sep=''))
    par(mar=c(4, 4, 4, 4))
    plot(x = x, y = dens, axes=F, xlim=c(x_min,x_max), ylim=c(0,max(dens,na.rm=TRUE)), 
         xlab="", ylab="", type="l", main="", col="gray90", lty=1, lwd = thickness)
    polygon(c(x, rev(x)), c(dens, rep(0,length(dens))), col = "gray90", lty = 0)
    par(new=T)
    plot(x = x, y = out_list[[plot]], type='l', axes=F, xlim=c(x_min,x_max), ylim=c(y_min,y_max), 
         xlab="", ylab="",col="blue", main="", lwd = thickness, cex.lab = font_size, cex.axis = font_size)
    grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
    if (plot == 'value') lines(x = x, y = res_val, col = 'red', lty = 2, lwd = thickness)
    if (plot == 'debt') lines(x = x, y = debt_commit, col = 'green', lty = 2, lwd = thickness)
    if (plot == 'spread') {
      lines(x = x, y = spread_commit, col = 'green', lty = 2, lwd = thickness)
      lines(x = x, y = risk_premium, col = 'red', lty = 5, lwd = thickness)
    }
    axis(2, col="black",lwd=2, cex.axis = axis_size)
    # mtext(2, text = legend[plot], line = 2, cex = font_size)
    title(ylab = legend[plot], line = 2.5, cex.lab = font_size)
    if (plot == 'cons_vol') abline(h = sigma, lty = 2, col = 'grey', lwd = thickness)
    if ((plot == 'ca')|(plot == 'g')|(plot == 'drift')) abline(h = 0, lty = 2, col = 'black', lwd = thickness)
    #======================== 
    # temporary additional line
    if (plot == 'cons') abline(h = c_post_default, lty = 2, col = 'grey', lwd = thickness)
    #======================== 
    abline(v = x_cutoff, lty = 3, col = 'purple', lwd = thickness)
    abline(v = theta*x_cutoff, lty = 4, col = 'orange', lwd = thickness)
    axis(1, col="black", lwd=2, cex.axis = axis_size)
    mtext(1,text=paste("debt-to-income x",sep=""),line=2.5,cex = font_size)
    box(which = "plot", bty = "l", lwd = thickness)
    dev.off()
  }
  
  # plot all functions of that will be in the main paper in a .eps format
  for (plot in c('debt','value','g','spread')){
    y_min <- min(out_list[[plot]], na.rm = TRUE)
    y_max <- max(out_list[[plot]], na.rm = TRUE)
    if (plot == 'g') {
      y_max <- min(1,y_max)
      y_min <- max(-1,y_min)
    }
    if (plot == 'spread') {
      y_min = 0
      y_max = min(.05,y_max)
    }
    if (plot == 'debt') {
      y_min = 0
      y_max = 1
    }
    
    setEPS()
    postscript(paste(plot,".eps",sep=''))
    par(mar=c(4, 4, 4, 4))
    plot(x = x, y = dens, axes=F, xlim=c(x_min,x_max), ylim=c(0,max(dens,na.rm=TRUE)), 
         xlab="", ylab="", type="l", main="", col="gray90", lty=1, lwd = thickness)
    polygon(c(x, rev(x)), c(dens, rep(0,length(dens))), col = "gray90", lty = 0)
    par(new=T)
    plot(x = x, y = out_list[[plot]], type='l', axes=F, xlim=c(x_min,x_max), ylim=c(y_min,y_max), 
         xlab="", ylab="",col="blue", main="", lwd = thickness, cex.lab = font_size, cex.axis = font_size)
    grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
    if (plot == 'value') lines(x = x, y = res_val, col = 'red', lty = 2, lwd = thickness)
    if (plot == 'debt') lines(x = x, y = debt_commit, col = 'green', lty = 2, lwd = thickness)
    if (plot == 'spread') {
      lines(x = x, y = spread_commit, col = 'green', lty = 2, lwd = thickness)
      lines(x = x, y = risk_premium, col = 'red', lty = 5, lwd = thickness)
    }
    axis(2, col="black",lwd=2, cex.axis = axis_size)
    title(ylab = legend[plot], line = 2.5, cex.lab = font_size)
    if (plot == 'cons_vol') abline(h = sigma, lty = 2, col = 'grey', lwd = thickness)
    if ((plot == 'ca')|(plot == 'g')|(plot == 'drift')) abline(h = 0, lty = 2, col = 'black', lwd = thickness)
    #======================== 
    # temporary additional line
    if (plot == 'cons') abline(h = c_post_default, lty = 2, col = 'grey', lwd = thickness)
    #======================== 
    abline(v = x_cutoff, lty = 3, col = 'purple', lwd = thickness)
    abline(v = theta*x_cutoff, lty = 4, col = 'orange', lwd = thickness)
    axis(1, col="black", lwd=2, cex.axis = axis_size)
    mtext(1,text=paste("debt-to-income x",sep=""),line=2.5,cex = font_size)
    box(which = "plot", bty = "l", lwd = thickness)
    dev.off()
  }
  
  # export data in data table
  out_list[['x']]  <- x
  out_df           <- as.data.frame(matrix(0, ncol = length(out_list)+3, nrow = n))
  names(out_df)    <- c(names(out_list),'v_def','debt_commit','spread_commit')
  for (dummy in names(out_df)){
    out_df[,dummy] <- out_list[[dummy]]
  }
  out_df[,'v_def']         <- res_val
  out_df[,'debt_commit']   <- debt_commit
  out_df[,'spread_commit'] <- spread_commit

  # export all policies and stationary distribution into csv file
  write.table(out_df, file = paste("plots.csv",sep=""), eol = "\n", na = "NA", fileEncoding = "",
              row.names = FALSE)
}

plot_sensi_fn <- function(param, plotting_param, sensi_param, sensi_vec){
  # function that plots all useful plots
  # this allows us to visualize sensitivity w.r.t. a given model parameter
  # sensi param is the name of the parameter to shock
  # sensi_vec is the range of values the parameter takes

  # plotting parameters
  n              <- plotting_param['n']
  thickness      <- plotting_param['thick']
  font_size      <- plotting_param['font_size']
  axis_size      <- plotting_param['axis_size']
  
  # compute number of sensi
  n_sensi <- length(sensi_vec)
  
  # declare cutoffs
  x_cutoff   <- array(NA,n_sensi)
  z_cutoff   <- array(NA,n_sensi)
  x_a        <- array(NA,n_sensi)
  x_0        <- array(NA,n_sensi)
  x_curr_act <- array(NA,n_sensi)
  x_reinject <- array(NA,n_sensi)
  z_reinject <- array(NA,n_sensi)
  
  # declare functions
  out_list   <- list()
  x          <- list()
  param_list <- list()

  # for each parameter set, compute relevant values
  for (i_sensi in (1:n_sensi)){
    # set parameter appropriately
    param[sensi_param]    <- sensi_vec[i_sensi]
    param_list[[i_sensi]] <- param
    
    # compute cutoff and other key cutoffs
    x_cutoff[i_sensi]   <- no_commit_cutoff_fn(param)
    z_cutoff[i_sensi]   <- (param['kappa']+param['m'])*x_cutoff[i_sensi]
    x_a[i_sensi]        <- no_commit_attraction_pt_fn(param)
    x_0[i_sensi]        <- no_commit_zero_drift_pt_fn(param)
    x_curr_act[i_sensi] <- no_commit_zero_ca_pt_fn(param)
    x_reinject[i_sensi] <- param['theta']*x_cutoff[i_sensi]
    z_reinject[i_sensi] <- (param['kappa']+param['m'])*x_reinject[i_sensi]
    
    # create grid
    x[[i_sensi]]        <- seq(from = .1, to = x_cutoff[i_sensi], length.out = n)
    
    # construct grid of output vectors
    out_list[[i_sensi]]                <- list()
    out_list[[i_sensi]][['debt']]      <- no_commit_debt_px_fn(param, x[[i_sensi]])
    out_list[[i_sensi]][['value']]     <- no_commit_value_fn(param, x[[i_sensi]])
    out_list[[i_sensi]][['g']]         <- no_commit_g_fn(param, x[[i_sensi]])
    out_list[[i_sensi]][['drift']]     <- no_commit_drift_fn(param, x[[i_sensi]])
    out_list[[i_sensi]][['cons']]      <- no_commit_c_fn(param, x[[i_sensi]])
    out_list[[i_sensi]][['ca']]        <- 1 - no_commit_c_fn(param, x[[i_sensi]])
    out_list[[i_sensi]][['spread']]    <- no_commit_credit_spread_fn(param, x[[i_sensi]])
    out_list[[i_sensi]][['cons_vol']]  <- no_commit_cons_vol_fn(param, x[[i_sensi]])
    dens_fn                            <- no_commit_ergodic_dist_fn(param)[['density']]
    out_list[[i_sensi]][['dens']]      <- dens_fn(x[[i_sensi]])
  }
  
  # name plots
  legend        <- c('debt price', 'value function', 'issuance policy', 'drift', 'consumption-to-income', 
                     'current account-to-income', 'credit spread', 'consumption volatility', 'density')
  names(legend) <- names(out_list[[1]])
  
  # color of lines
  legend_col  <- colours()[10+(1:n_sensi)*5]
  legend_name <- paste0(sensi_param,' = ', sensi_vec)

  # plot all functions of interest
  x_max <- max(x_cutoff)
  x_min <- .5*min(x_reinject)
  z_max <- max(z_cutoff)
  z_min <- .5*min(z_reinject)
  for (plot in names(out_list[[1]])){
    y_min <- min(out_list[[1]][[plot]], na.rm = TRUE)
    y_max <- max(out_list[[1]][[plot]], na.rm = TRUE)
    for (i_sensi in (2:n_sensi)){
      y_min <- min(y_min, out_list[[i_sensi]][[plot]], na.rm = TRUE)
      y_max <- max(y_max, out_list[[i_sensi]][[plot]], na.rm = TRUE)
    }
    if ((plot == 'ca')|(plot == 'g')|(plot == 'drift')) {
      y_max <- 1
      y_min <- -1
    }
    if (plot == 'cons')     y_max = 2
    if (plot == 'cons_vol') y_min = 0
    if (plot == 'spread') {
      y_min = 0
      y_max = .05
    }
    if (plot == 'debt') {
      y_min = 0
      y_max = 1
    }
    # plot as a function of x
    png(paste0("comp_stat_",sensi_param,"_",plot,"_vs_x.png"))
    par(mar=c(4, 4, 4, 4))
    plot(x = x[[1]], y = out_list[[1]][[plot]], type='l', axes=F, xlim=c(x_min,x_max), ylim=c(y_min,y_max), lty = 1, 
         xlab="", ylab="",col=legend_col[1], main="", lwd = thickness, cex.lab = font_size, cex.axis = font_size)
    grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
    for (i_sensi in (2:n_sensi)) lines(x = x[[i_sensi]], y = out_list[[i_sensi]][[plot]], col=legend_col[i_sensi], lwd=thickness, lty=i_sensi)
    axis(2, col="black",lwd=2, cex.axis = axis_size)
    mtext(2,text = legend[plot], line = 2, cex = font_size)
    if ((plot == 'ca')|(plot == 'g')|(plot == 'drift')) abline(h = 0, lty = 2, col = 'grey', lwd = thickness)
    axis(1, col="black", lwd=2, cex.axis = axis_size)
    mtext(1,text=paste("debt-to-income x",sep=""),line=2.5,cex = font_size)
    legend("topright", col = legend_col, bty="n", lty = 1:n_sensi, cex = font_size, legend = legend_name, lwd = thickness)
    box(which = "plot", bty = "l", lwd = thickness)
    dev.off()
    
    # plot as a function of z
    png(paste0("comp_stat_",sensi_param,"_",plot,"_vs_z.png"))
    par(mar=c(4, 4, 4, 4))
    plot(x = (param_list[[1]]['m'] + param_list[[1]]['kappa'])*x[[1]], y = out_list[[1]][[plot]], type='l', axes=F, xlim=c(z_min,z_max), ylim=c(y_min,y_max), lty = 1, 
         xlab="", ylab="",col=legend_col[1], main="", lwd = thickness, cex.lab = font_size, cex.axis = font_size)
    grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
    for (i_sensi in (2:n_sensi)) lines(x = (param_list[[i_sensi]]['m'] + param_list[[i_sensi]]['kappa'])*x[[i_sensi]], y = out_list[[i_sensi]][[plot]], col=legend_col[i_sensi], lwd=thickness, lty=i_sensi)
    axis(2, col="black",lwd=2, cex.axis = axis_size)
    mtext(2,text = legend[plot], line = 2, cex = font_size)
    if ((plot == 'ca')|(plot == 'g')|(plot == 'drift')) abline(h = 0, lty = 2, col = 'grey', lwd = thickness)
    axis(1, col="black", lwd=2, cex.axis = axis_size)
    mtext(1,text=paste("debt coverage ratio z",sep=""),line=2.5,cex = font_size)
    legend("topright", col = legend_col, bty="n", lty = 1:n_sensi, cex = font_size, legend = legend_name, lwd = thickness)
    box(which = "plot", bty = "l", lwd = thickness)
    dev.off()
  }
}

sim_fn <- function(param_vec, sim_param, plotting_param){
  # function that computes random path for all variables of interest
  set.seed(908)
  
  # set max issuance rate
  max_drift <- 1000
  
  # number of years of simulation
  T_sim    <- sim_param['T']
  dt       <- sim_param['dt']
  Y_ini    <- sim_param['Y0']
  F_ini    <- sim_param['F0']
  
  # plotting parameters
  thickness <- plotting_param['thick']
  font_size <- plotting_param['font_size']
  axis_size <- plotting_param['axis_size']
  
  # download parameter theta
  alpha    <- param_vec['alpha']
  theta    <- param_vec['theta']
  mu       <- param_vec['mu']
  sigma    <- param_vec['sigma']
  m        <- param_vec['m']
  
  # compute time vector
  t_vec    <- seq(from = 0, to = T_sim, by = dt)
  dt_vec   <- diff(t_vec)
  n_t      <- length(dt_vec)
  x_ini    <- F_ini/Y_ini
  y_ini    <- Y_ini/F_ini
  
  # output vector
  Y_vec     <- array(Y_ini,n_t)
  Y_vec2    <- array(Y_ini,n_t)
  F_vec     <- array(F_ini,n_t)
  F_vec2    <- array(F_ini,n_t)
  x_vec     <- array(x_ini,n_t)
  y_vec     <- array(y_ini,n_t)
  def_vec   <- array(0,n_t)
  def_vec2  <- array(0,n_t)
  drift_x   <- array(0,n_t)
  drift_y   <- array(0,n_t)
  
  # compute cutoff and attraction point
  x_bar     <- no_commit_cutoff_fn(param_vec)
  y_bar     <- 1 / x_bar
  att_pt_x  <- no_commit_attraction_pt_fn(param_vec)
  att_pt_y  <- no_commit_attraction_pt_fn2(param_vec)

  # simulate shocks
  shock_vec <- rnorm(n_t)
  for (t in (1:(n_t-1))){
    drift_x[t] <- min(max_drift,no_commit_drift_fn(param_vec, x_vec[t]))
    drift_y[t] <- min(max_drift,no_commit_drift_fn2(param_vec, y_vec[t]))
    diff_x       <- - sigma * x_vec[t] * shock_vec[t]
    diff_y       <- sigma * y_vec[t] * shock_vec[t]
    x_vec[t+1]   <- x_vec[t] + drift_x[t]*dt + diff_x*sqrt(dt)
    y_vec[t+1]   <- y_vec[t] + drift_y[t]*dt + diff_y*sqrt(dt)
    def_vec[t+1] <- 0*(x_vec[t+1]<x_bar) + 1*(x_vec[t+1]>=x_bar)
    def_vec2[t+1]<- 0*(y_vec[t+1]>y_bar) + 1*(y_vec[t+1]<=y_bar)
    x_vec[t+1]   <- x_vec[t+1]*(x_vec[t+1]<x_bar) + theta * x_bar *(x_vec[t+1]>=x_bar)
    y_vec[t+1]   <- y_vec[t+1]*(y_vec[t+1]>y_bar) + (1/theta) * y_bar *(y_vec[t+1]<=y_bar)
    Y_vec[t+1]   <- Y_vec[t]*(1 + mu * dt + sigma * shock_vec[t]*sqrt(dt)) * (1 + (alpha - 1)*def_vec[t+1])
    Y_vec2[t+1]  <- Y_vec2[t]*(1 + mu * dt + sigma * shock_vec[t]*sqrt(dt)) * (1 + (alpha - 1)*def_vec2[t+1])
    F_vec[t+1]   <- Y_vec[t+1]*x_vec[t+1]
    F_vec2[t+1]  <- Y_vec2[t+1]/y_vec[t+1]
  }
  
  # plot time series of debt-to-income ratio
  if (is.nan(att_pt_x))  {
    y_min <- .9*min(min(x_vec),theta*x_bar)
  } else {
    y_min <- .9*min(theta*x_bar,att_pt_x,min(x_vec))
  }
  y_max <- 1.1*x_bar
  
  # create pdf
  pdf("debt_to_income.pdf")
  plot(x = t_vec[-1], y = x_vec, lwd = thickness, ylim = c(y_min,y_max),
       bty = 'n', type = 'l', ylab = 'debt-to-income x', xlab = 'time', col = 'blue', cex.lab = font_size, cex.axis = font_size)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  plot_coord <- par('usr')
  abline(h = x_bar, col = 'red', lty = 2, lwd = thickness)
  abline(h = att_pt_x, col = 'green', lty = 3, lwd = thickness)
  # abline(h = theta*x_bar, col = 'black', lty = 4, lwd = thickness)
  text(y = x_bar + .03*(plot_coord[4]-plot_coord[3]), 
       x = plot_coord[1] + .05*(plot_coord[2]-plot_coord[1]), 
       expression(bar(x)), cex = font_size, col = 'black') 
  # text(y = theta*x_bar - .03*(plot_coord[4]-plot_coord[3]), 
  #      x = plot_coord[1] + .05*(plot_coord[2]-plot_coord[1]), 
  #      expression(paste(theta,bar(x))), cex = font_size, col = 'black') 
  text(y = att_pt_x + .03*(plot_coord[4]-plot_coord[3]), 
       x = plot_coord[1] + .05*(plot_coord[2]-plot_coord[1]), 
       expression(x[a]), cex = font_size, col = 'black') 
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()

  # create eps
  setEPS()
  postscript("debt_to_income.eps")
  plot(x = t_vec[-1], y = x_vec, lwd = thickness, ylim = c(y_min,y_max),
       bty = 'n', type = 'l', ylab = 'debt-to-income x', xlab = 'time', col = 'blue', cex.lab = font_size, cex.axis = font_size)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  plot_coord <- par('usr')
  abline(h = x_bar, col = 'red', lty = 2, lwd = thickness)
  abline(h = att_pt_x, col = 'green', lty = 3, lwd = thickness)
  # abline(h = theta*x_bar, col = 'black', lty = 4, lwd = thickness)
  text(y = x_bar + .03*(plot_coord[4]-plot_coord[3]), 
       x = plot_coord[1] + .05*(plot_coord[2]-plot_coord[1]), 
       expression(bar(x)), cex = font_size, col = 'black') 
  # text(y = theta*x_bar - .03*(plot_coord[4]-plot_coord[3]), 
  #      x = plot_coord[1] + .05*(plot_coord[2]-plot_coord[1]), 
  #      expression(paste(theta,bar(x))), cex = font_size, col = 'black') 
  text(y = att_pt_x + .03*(plot_coord[4]-plot_coord[3]), 
       x = plot_coord[1] + .05*(plot_coord[2]-plot_coord[1]), 
       expression(x[a]), cex = font_size, col = 'black') 
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # plot time series of debt-to-income ratio
  y_min <- .9*y_bar
  y_max <- 1.1*max(theta*y_bar, att_pt_y, quantile(y_vec,.99))
  pdf("income_to_debt.pdf")
  plot(x = t_vec[-1], y = y_vec, lwd = thickness, ylim = c(y_min, y_max),
       bty = 'n', type = 'l', ylab = 'income-to-debt y', xlab = 'time', col = 'blue', cex.lab = font_size, cex.axis = font_size)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  plot_coord <- par('usr')
  abline(h = y_bar, col = 'red', lty = 2, lwd = thickness)
  abline(h = att_pt_y, col = 'green', lty = 3, lwd = thickness)
  # abline(h = y_bar/theta, col = 'blue', lty = 4, lwd = thickness)
  text(y = y_bar + .03*(plot_coord[4]-plot_coord[3]), 
       x = plot_coord[1] + .05*(plot_coord[2]-plot_coord[1]), 
       expression(bar(y)), cex = font_size, col = 'black') 
  # text(y = y_bar/theta + .03*(plot_coord[4]-plot_coord[3]), 
  #      x = plot_coord[1] + .05*(plot_coord[2]-plot_coord[1]), 
  #      expression(bar(y)/theta), cex = font_size, col = 'black') 
  text(y = att_pt_y + .03*(plot_coord[4]-plot_coord[3]), 
       x = plot_coord[1] + .05*(plot_coord[2]-plot_coord[1]), 
       expression(y[a]), cex = font_size, col = 'black') 
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # plot time series of debt and income
  pdf("debt_and_income.pdf")
  plot(x = t_vec[-1], y = F_vec, lwd = thickness, ylim = c(0,1.1*max(F_vec,Y_vec)), lty = 2, 
       bty = 'n', type = 'l', ylab = 'debt notional F and income Y', xlab = 'time', col = 'blue', 
       cex.lab = font_size, cex.axis = font_size)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  plot_coord <- par('usr')
  lines(x = t_vec[-1], y = Y_vec, type = 'l', col = 'red', lty = 1,lwd = thickness)
  legend("bottomright", legend = c('F(t)','Y(t)'), 
         col = c('blue', 'red'), bty="n", lty = c(2,1), cex = font_size,lwd = thickness)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # plot time series of debt and income in .eps format
  setEPS()
  postscript("debt_and_income.eps")
  plot(x = t_vec[-1], y = F_vec, lwd = thickness, ylim = c(0,1.1*max(F_vec,Y_vec)), lty = 2, 
       bty = 'n', type = 'l', ylab = 'debt notional F and income Y', xlab = 'time', col = 'blue', 
       cex.lab = font_size, cex.axis = font_size)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  plot_coord <- par('usr')
  lines(x = t_vec[-1], y = Y_vec, type = 'l', col = 'red', lty = 1,lwd = thickness)
  legend("bottomright", legend = c('F(t)','Y(t)'), 
         col = c('blue', 'red'), bty="n", lty = c(2,1), cex = font_size,lwd = thickness)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # export data in data table
  out_df           <- as.data.frame(matrix(0, ncol = 4, nrow = n_t))
  names(out_df)    <- c('t','x','F','Y')
  out_df[,'t']     <- t_vec[-1] 
  out_df[,'x']     <- x_vec 
  out_df[,'Y']     <- Y_vec
  out_df[,'F']     <- F_vec

  # export simulation into csv file
  write.table(out_df, file = paste("simulation.csv",sep=""), eol = "\n", na = "NA", fileEncoding = "", 
              row.names = FALSE)
}

#====================================================================
# Functions for model where sovereign commits to G(t) = gamma F(t) dt
#====================================================================

commit_debt_root_fn <- function(param_vec){
  # compute root xi for debt function with commitment
  # assumption: the government can commit to issue gamma F(t) dt per unit of time
  delta   <- param_vec["delta"]
  m       <- param_vec["m"]
  mu      <- param_vec["mu"]
  sigma   <- param_vec["sigma"]
  gamma   <- param_vec["gamma"]
  r       <- param_vec["r"]
  nu      <- param_vec["nu"]
  if (is.na(gamma)) gamma <- 0 # in case we forgot to specify gamma, the issuance rate
  
  # compute root of characteristic polynomial for debt HJB
  a       <- 0.5*sigma^2
  b       <- - (m + mu - nu * sigma - gamma - 0.5*sigma^2)
  c       <- - (r + m)
  discrim <- b^2 - 4*a*c
  xi      <- (- b + sqrt(discrim))/(2*a)
  
  return(xi)
}

commit_debt_root2_fn <- function(param){
  # compute root xi for debt function with commitment -- generate the 2 roots
  # assumption: the government can commit to issue gamma F(t) dt per unit of time
  delta    <- param["delta"]
  m        <- param["m"]
  mu       <- param["mu"]
  sigma    <- param["sigma"]
  gamma    <- param["gamma"]
  r        <- param["r"]
  if (is.na(gamma)) gamma <- 0 # in case we forgot to specify gamma, the issuance rate
  
  # compute roots
  a          <- 0.5*sigma^2
  b          <- - (m + mu - gamma - 0.5*sigma^2)
  c          <- - (r + m)
  discrim    <- b^2 - 4*a*c
  eta_2      <- (- b + sqrt(discrim))/(2*a)
  eta_1      <- (- b - sqrt(discrim))/(2*a)
  eta        <- c(eta_1,eta_2)
  names(eta) <- c("eta_1", "eta_2")
  return(eta)
}

commit_value_root_fn <- function(param_vec){
  # compute root xi for value function with commitment
  # assumption: the government can commit to issue gamma F(t) dt per unit of time
  delta   <- param_vec["delta"]
  m       <- param_vec["m"]
  mu      <- param_vec["mu"]
  sigma   <- param_vec["sigma"]
  gamma   <- param_vec["gamma"]
  if (is.na(gamma)) gamma <- 0 # in case we forgot to specify gamma, the issuance rate
  
  # compute root of characteristic polynomial for value function HJB
  a       <- 0.5*sigma^2
  b       <- - (m + mu - gamma + 0.5*sigma^2)
  c       <- - (delta - mu)
  discrim <- b^2 - 4*a*c
  xi      <- (- b + sqrt(discrim))/(2*a)
  
  return(xi)
}

commit_cutoff_fn <- function(param_vec){
  # compute optimal default cutoff for case with commitment
  # assumption: the government can commit to issue gamma F(t) dt per unit of time
  r       <- param_vec["r"]
  delta   <- param_vec["delta"]
  kappa   <- param_vec["kappa"]
  m       <- param_vec["m"]
  mu      <- param_vec["mu"]
  sigma   <- param_vec["sigma"]
  alpha   <- param_vec["alpha"]
  theta   <- param_vec["theta"]
  gamma   <- param_vec["gamma"]
  if (is.na(gamma)) gamma <- 0 # in case we forgot to specify gamma, the issuance rate
  
  # compute root of characteristic polynomial for debt
  xi_d    <- commit_debt_root_fn(param_vec)
  
  # compute root of characteristic polynomial for value function
  xi_v    <- commit_value_root_fn(param_vec)
  
  # compute default boundary
  x_bar   <- xi_v/(((r + m - gamma)/(delta + m - gamma) - gamma/(r + gamma - delta))*(xi_v-1) + 
                     (gamma/(r + gamma - delta))*xi_d) * 
    ((r + m)/(kappa + m)) * ((1 - alpha)/(1 - alpha * theta)) / (delta - mu)
  
  return(x_bar)
}

commit_debt_px_fn <- function(param_vec, x){
  # compute debt price for case with commitment
  # assumption: the government can commit to issue gamma F(t) dt per unit of time
  delta   <- param_vec["delta"]
  kappa   <- param_vec["kappa"]
  m       <- param_vec["m"]
  mu      <- param_vec["mu"]
  sigma   <- param_vec["sigma"]
  alpha   <- param_vec["alpha"]
  theta   <- param_vec["theta"]
  gamma   <- param_vec["gamma"]
  r       <- param_vec["r"]
  nu      <- param_vec["nu"]
  if (is.na(gamma)) gamma <- 0 # in case we forgot to specify gamma, the issuance rate
  
  # compute root of characteristic polynomial for debt
  xi_d    <- commit_debt_root_fn(param_vec)
  
  # compute default boundary
  x_bar   <- commit_cutoff_fn(param_vec)
  
  # compute debt price
  debt    <- (kappa + m)/(r + m) * 
    (1 - (1 - alpha*theta)/(1 - alpha*(theta^(xi_d+1))) * (x/x_bar)^xi_d)
  
  # take care of situations where x > x_bar -- compute number if consecutive restructurings
  n       <- pmax(0,1 + floor((log(x_bar) - log(x))/log(theta)))
  
  # compute debt price in case x > x_bar
  debt_def<- ((alpha*theta)^n) * ((kappa + m)/(r + m) * 
                                    (1 - (1 - alpha*theta)/(1 - alpha*(theta^(xi_d+1))) * ((theta^n)*x/x_bar)^xi_d))
  # return debt price
  out <- debt*(n<=0) + debt_def*(n>0)
  
  return(out)
}

commit_value_fn <- function(param_vec, x){
  # compute value function for case with commitment
  # assumption: the government can commit to issue gamma F(t) dt per unit of time
  delta   <- param_vec["delta"]
  kappa   <- param_vec["kappa"]
  m       <- param_vec["m"]
  mu      <- param_vec["mu"]
  sigma   <- param_vec["sigma"]
  alpha   <- param_vec["alpha"]
  theta   <- param_vec["theta"]
  gamma   <- param_vec["gamma"]
  r       <- param_vec["r"]
  if (is.na(gamma)) gamma <- 0 # in case we forgot to specify gamma, the issuance rate
  
  # compute root of characteristic polynomial for debt
  xi_d    <- commit_debt_root_fn(param_vec)
  # compute root of characteristic polynomial for value function
  xi_v    <- commit_value_root_fn(param_vec)
  # compute default boundary
  x_bar   <- commit_cutoff_fn(param_vec)
  # compute value function
  value   <- 1 / (delta - mu) * 
    (1 - (1 - alpha)/(1 - alpha*(theta^xi_v)) * (x/x_bar)^(xi_v)) - 
    x * (kappa + m)*(r + m - gamma)/(delta + m - gamma)/(r + m) * 
    (1 - (1 - alpha*theta)/(1 - alpha*(theta^xi_v)) * (x/x_bar)^(xi_v - 1))+
    x * (kappa + m)*gamma/(r + gamma - delta)/(r + m) * 
    ((1 - alpha*theta)/(1 - alpha*(theta^(xi_d+1))) * (x/x_bar)^(xi_d) - 
       (1 - alpha*theta)/(1 - alpha*(theta^xi_v)) * (x/x_bar)^(xi_v - 1))
  
  return(value)
}

commit_credit_spread_fn <- function(param_vec, x){
  # compute credit spread function for case with commitment
  # assumption: the government can commit to issue gamma F(t) dt per unit of time
  delta   <- param_vec["delta"]
  kappa   <- param_vec["kappa"]
  m       <- param_vec["m"]
  mu      <- param_vec["mu"]
  sigma   <- param_vec["sigma"]
  alpha   <- param_vec["alpha"]
  theta   <- param_vec["theta"]
  r       <- param_vec["r"]
  nu      <- param_vec["nu"]
  gamma   <- param_vec["gamma"]
  if (is.na(gamma)) gamma <- 0 # in case we forgot to specify gamma, the issuance rate
  
  # compute root of characteristic polynomial for debt
  xi_d    <- commit_debt_root_fn(param_vec)
  # compute default boundary
  x_bar   <- commit_cutoff_fn(param_vec)
  # compute credit spreads  
  spread    <- (m + r) / (1 - (1 - alpha*theta)/(1 - alpha*theta^(xi_d+1)) * (x/x_bar)^(xi_d)) - (r + m)
  
  return(spread)
}

#=====================================================================
# Citizens' Welfare Function -- Citizens more patients than government
#=====================================================================

no_commit_hat_value_fn <- function(param, num_param){
  # compute value function for citizens in case without commitment
  # compute for given assumed boundaries

  # load numerical parameters
  n_x       <- num_param['n_x']
  dt        <- num_param['dt']
  max_iter  <- num_param['max_iter']
  tol       <- num_param['tol']

  # load parameters
  delta     <- param["delta"]
  hat_delta <- param["hat_delta"]
  kappa     <- param["kappa"]
  m         <- param["m"]
  mu        <- param["mu"]
  sigma     <- param["sigma"]
  alpha     <- param["alpha"]
  theta     <- param["theta"]
  r         <- param["r"]
  
  # compute default boundary without commitment
  x_bar    <- no_commit_cutoff_fn(param)
  
  # set x_min so that we do not need to deal with singularity
  x_min    <- .01

  # create grid
  x     <- seq(from = x_min, to = x_bar, length.out = n_x)
  dx    <- x[2]-x[1]
  
  # compute issuance and consumption policies, as well as original guess for v_n
  g      <- no_commit_g_fn(param, x)
  cons   <- no_commit_c_fn(param, x)
  orig_v <- no_commit_value_fn(param, x)
  v_n    <- orig_v
  
  # compute drift rate (invariant across iterations)
  drift  <- g - (mu + m)*x
  
  # construct matrix A
  A      <- Matrix(0, nrow = n_x, ncol = n_x, sparse = TRUE)
  
  # diagonal elements for A
  locs     <- cbind(seq.int(n_x),seq.int(n_x))
  A[locs]  <- -(abs(drift)/dx + (sigma*x/dx)^2)
  
  # lower diagonal elements -- probability to go to lower x
  locs     <- cbind((2:n_x),(1:(n_x-1)))
  A[locs]  <- -pmin(0,drift[-1])/dx + .5*(sigma*x[-1]/dx)^2
  
  # upper diagonal elements -- probability to go to higher x
  locs     <- cbind((1:(n_x-1)),(2:n_x))
  A[locs]  <- pmax(0,drift[-n_x])/dx + .5*(sigma*x[-n_x]/dx)^2
  
  # impose reflection at x = 0 and x = x_bar
  A[1,1]     <- -sum(A[1,-1])
  A[n_x,n_x] <- -sum(A[n_x,-n_x])
  
  # construct B matrix
  B              <- (1/dt + hat_delta - mu)*Diagonal(n_x) - A
  
  # boundary condition at x = x_bar
  idx_reentry        <- which.min(abs(x-theta*x_bar))
  B[n_x,]            <- 0
  B[n_x,n_x]         <- 1
  B[n_x,idx_reentry] <- -alpha
  
  # constant part of the right hand side + take care of boundary conditions
  RHS_cst            <- cons

  # initialization
  err     <- 1      # initialize error for value function
  err_vec <- c(1)
  iter    <- 1      # initialize iteration #
  # compute citizens' welfare value
  while ((err>tol)&(iter<max_iter)){
    # if error greater than tolerance
    # assign most recent value to v
    v        <- v_n
    
    # enter right hand side, including correction for boundary conditions
    RHS      <- RHS_cst + v/dt
    RHS[n_x] <- 0
    
    # solve linear system
    v_n  <- c(array(solve(B,RHS)))

    # compute error
    err      <- max(abs(v_n-v)/dt)
    err_vec  <- c(err_vec,err)
    
    # update iteration #
    iter     <- iter + 1
  }
  
  # interpolate at zero
  hat_v_0  <- v_n[1] - x_min*(v_n[2]-v_n[1])/dx
  x        <- c(0,x)
  v        <- c(hat_v_0,v_n)
  
  # smooth function
  hat_v_fn <- approxfun(x, v, method = "linear")
  
  # output
  output   <- list(hat_v_fn,hat_v_0,iter,err)
  names(output) <- c("hat_v", "hat_v0", "n_iter", "error")
  
  return(output)
}

#===================================================
# Value function Poisson opportunities to sell bonds
# no commitment
#===================================================

Poisson_no_commit_fn <- function(param, num_param, plotting_param){
  # compute value function for case where financing opportunities arrive at Poisson times 

  # setup for verification purposes
  dump_folder<- "C:/Users/ft.fi/Dropbox/Research/Sovereign Defaults and Commitment/Code/Dump"

  # discretization points for numerical solution
  n          <- num_param['n_x']        # number of discretization points
  dt         <- num_param['dt']         # time interval for false transient method
  max_iter   <- num_param['max_iter']   # max number of iterations
  tol        <- num_param['tol']        # tolerance
  print_flag <- num_param['print_flag'] # print flag for error checking purposes
  
  # load parameters
  delta     <- param["delta"]
  lambda    <- param["lambda"]
  kappa     <- param["kappa"]
  m         <- param["m"]
  mu        <- param["mu"]
  nu        <- param["nu"]
  sigma     <- param["sigma"]
  alpha     <- param["alpha"]
  theta     <- param["theta"]
  r         <- param["r"]
  
  # plotting parameters
  thickness      <- plotting_param['thick']
  font_size      <- plotting_param['font_size']
  axis_size      <- plotting_param['axis_size']
  
  # compute default boundary without commitment -- this is our initial guess default boundary
  x_bar_nc  <- no_commit_cutoff_fn(param)
  x_bar_c   <- commit_cutoff_fn(param)
  x_bar     <- x_bar_nc

  # set x_min so that we do not need to deal with singularity
  x_min    <- .01
  
  # compute debt price between x = 0 and x = 2*x_star
  x        <- seq(from = x_min, to = 1.5*x_bar, length.out = n)
  dx       <- x[2]-x[1]
  dx2      <- dx^2
  i_bar    <- which.min(abs(x - x_bar)) 
  x_bar    <- x[i_bar]

  # compute debt price assuming commitment to never issue debt in future
  d_c     <- commit_debt_px_fn(param, x)
  
  # compute initial guess for debt price, value function and issuance policy
  # initial guess rely on no-commitment value function
  v_nc    <- no_commit_value_fn(param, x)
  d_nc    <- no_commit_debt_px_fn(param, x)
  g_nc    <- no_commit_g_fn(param, x)
  v_n     <- v_nc
  d_n     <- d_nc
  jump_to <- x
  i_next  <- apply(as.matrix(1:n), 1, FUN = function(i) return(which.min(abs(x - jump_to[i]))))
  jump_val<- 0

  # compute drift rate (invariant across iterations)
  drift_v <- - (mu + m)*x
  drift_d <- - (mu + m - sigma^2 - nu*sigma)*x
  
  # construct matrix A_v and A_d
  A_v   <- Matrix(0, nrow = n, ncol = n, sparse = TRUE)
  A_d   <- Matrix(0, nrow = n, ncol = n, sparse = TRUE)
  
  # diagonal elements for A
  locs     <- cbind(seq.int(n),seq.int(n))
  A_v[locs]<- -(abs(drift_v)/dx + (sigma*x)^2/dx2)
  A_d[locs]<- -(abs(drift_d)/dx + (sigma*x)^2/dx2)
  
  # lower diagonal elements -- probability to go to lower x
  locs     <- cbind((2:n),(1:(n-1)))
  A_v[locs]<- -pmin(0,drift_v[-1])/dx + .5*(sigma*x[-1])^2/dx2
  A_d[locs]<- -pmin(0,drift_d[-1])/dx + .5*(sigma*x[-1])^2/dx2
  
  # upper diagonal elements -- probability to go to higher x
  locs     <- cbind((1:(n-1)),(2:n))
  A_v[locs]<- pmax(0,drift_v[-n])/dx + .5*(sigma*x[-n])^2/dx2
  A_d[locs]<- pmax(0,drift_d[-n])/dx + .5*(sigma*x[-n])^2/dx2
  
  # impose reflection at boundaries
  A_v[1,1] <- -sum(A_v[1,-1])
  A_v[n,n] <- -sum(A_v[n,-n])
  A_d[1,1] <- -sum(A_d[1,-1])
  A_d[n,n] <- -sum(A_d[n,-n])
  
  # construct B_v matrix
  B_v      <- (1/dt + delta - mu)*Diagonal(n) - A_v
  B_d      <- (1/dt + r + m)*Diagonal(n) - A_d
  
  # initialization
  err     <- 1      # initialize error for value function
  err_t   <- c(1)   # keep track of errors
  iter    <- 1      # initialize iteration #
  
  # start algo
  while ((err>tol)&(iter<max_iter)){
    # if error greater than tolerance
    v    <- v_n
    d    <- d_n
    
    # compute default option value for equity, and recovery value for debt at default
    v_interp <- splinefun(x, v, method = "natural")
    d_interp <- splinefun(x, d, method = "natural")
    v_def    <- alpha * v_interp(theta * x)
    d_def    <- alpha * theta * d_interp(theta * x)
    
    # right hand side for debt and value function
    payoff   <- 1 - (kappa+m)*x + jump_val
    RHS_d    <- kappa + m + lambda * (d[i_next] - d) + d / dt
    
    # impose value matching at x > x_bar for d
    RHS_d[i_bar:n] <- d_def[i_bar:n]
    
    # solve for debt price, imposing value matching at default boundary
    C_d                      <- B_d
    C_d[(i_bar:n),]          <- Matrix(0,nrow = (n-i_bar+1), ncol = n)
    C_d[(i_bar:n),(i_bar:n)] <- Diagonal(n-i_bar+1)
    
    # compute d(n)
    d_n   <- c(array(solve(C_d,RHS_d)))

    # frame as an LCP problem -- iteration n+1 as a function of iteration n
    # z(n+1) := v(n+1) - def_opt(n)   >= 0 --- always the option to default
    # B(n) z(n+1) + q(n)              >= 0 --- HJB in continuation region
    # q(n) = - flow payoff(n) - e(n) / dt + B(n) def_opt(n)
    # z(n+1)'(B(n)z(n+1) + q(n))  = 0 and minimize quadratic form
    q     <- - payoff - v/dt + c(array(B_v %*% v_def))
    
    # compute z(n+1) and thus e(n+1) using LCP solver
    z_guess   <- v - v_def
    lower_bd  <- rep(0,n)
    upper_bd  <- rep(Inf,n)
    z         <- LCP_fn(B_v,q,lower_bd,upper_bd,z_guess,0)
    v_n       <- z + v_def

    # keep track of LCP algo error
    LCP_error     <- abs(z*c(array(B_v %*% z + q)))

    # update x_high, and make sure x_bar is on the grid
    def_opt_id  <- which(v_n <= v_def + tol)
    if (length(def_opt_id) != 0){
      # update x_bar only if LCP algo did find area of state space with default exercise
      # otherwise don't update x_bar
      i_bar    <- min(def_opt_id)
      x_bar    <- x[i_bar]
    }
    
    # update jump policy
    i_next[1:i_bar]     <- apply(as.matrix(1:i_bar), 1, 
                                 FUN = function(i) return(which.max(v_n[1:i_bar] + (x[1:i_bar]-x[i])*d_n[1:i_bar])))
    i_next[(i_bar+1):n] <- (i_bar+1):n
    jump_to             <- x[i_next]
    jump_val            <- lambda * (v_n[i_next] + (x[i_next] - x) * d_n[i_next] - v_n)
    
    if (print_flag) {
      # plot for verification purposes just in case
      setwd(dump_folder)
      
      # create list to be plotted
      plots                <- list()
      plots[['value']]     <- v_n
      plots[['jump_val']]  <- jump_val
      plots[['debt']]      <- d_n
      plots[['jump_to']]   <- jump_to
      plots[['exp_iss']]   <- (jump_to - x)*lambda
      plots[['LCP_err']]   <- LCP_error

      # plot all functions of interest
      for (plt in names(plots)){
        y_min <- min(plots[[plt]])
        y_max <- max(plots[[plt]])
        png(paste(plt,"_iter_",iter,".png",sep=""))
        plot(x = x, y = plots[[plt]], bty = 'n', type = 'l', ylim = c(y_min,y_max), col = 'blue', xlab = 'debt-to-income', ylab = plt)
        if (plt == 'jump_to')   lines(x = x, y = x, lty = 2, col = 'green')
        if (plt == 'exp_iss')   lines(x = x, y = g_nc, lty = 2, col = 'green')
        if (plt == 'debt'){
          lines(x = x, y = d_nc, lty = 2, col = 'green')
          lines(x = x, y = d_c, lty = 3, col = 'purple')
        }
        if (plt == 'value')     lines(x = x, y = v_nc, lty = 2, col = 'green')
        abline(v = x_bar, col = 'red', lty = 2)
        abline(v = x_bar_nc, col = 'pink', lty = 3)
        dev.off()
      }
    }
    
    # compute error
    err      <- max(abs(v_n-v)/dt) 
    err_t    <- c(err_t,err)
    
    # update iteration #
    iter     <- iter + 1
  }
  
  # print number of iterations needed
  print(paste0('iterations needed: ',iter))
  
  # interpolate at zero
  v_0       <- v_n[1] - x_min*(v_n[2]-v_n[1])/dx
  d_0       <- d_n[1] - x_min*(d_n[2]-d_n[1])/dx
  jump_to_0 <- jump_to[1] - x_min*(jump_to[2]-jump_to[1])/dx
  jump_val_0<- jump_val[1] - x_min*(jump_val[2]-jump_val[1])/dx
  x         <- c(0,x)
  v         <- c(v_0,v_n)
  d         <- c(d_0,d_n)
  jump_to   <- c(jump_to_0,jump_to)
  jump_val  <- c(jump_val_0,jump_val)
  exp_iss   <- (jump_to - x)*lambda

  # find attraction point
  x_a       <- x[which.min(abs(exp_iss - m*x))]
  
  # generate plot for value function
  y_min <- min(v,v_nc)
  y_max <- max(v,v_nc)
  pdf('poisson_v.pdf')
  par(mar=c(4, 4, 4, 4))
  plot(x = x, y = v, type='l', axes=F, ylim=c(y_min,y_max), 
       xlab="", ylab="", col="blue", main="", lwd = thickness, cex.lab = font_size, cex.axis = font_size)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  lines(x = x[-1], y = v_nc, col = 'red', lty = 2, lwd = thickness)
  axis(2, col="black",lwd=2, cex.axis = axis_size)
  mtext(2,text = 'value function', line = 2, cex = font_size)
  axis(1, col="black", lwd=2, cex.axis = axis_size)
  mtext(1,text=paste("debt-to-income x",sep=""), line=2.5, cex = font_size)
  abline(v = x_bar, col = 'blue', lty = 4, lwd = thickness)
  abline(v = x_bar_nc, col = 'red', lty = 5, lwd = thickness)
  legend("topright", legend = c('poisson','no-commitment'), 
         col = c('blue', 'red'), bty="n", lty = c(1,2), cex = font_size, lwd = thickness)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # generate plot for value function differential between autarky and trade
  y_min <- min(v[-1]-v_nc,0)
  y_max <- max(v[-1]-v_nc,0)
  pdf('poisson_v_diff.pdf')
  par(mar=c(4, 4, 4, 4))
  plot(x = x[-1], y = v[-1]-v_nc, type='l', axes=F, ylim=c(y_min,y_max), 
       xlab="", ylab="", col="blue", main="", lwd = thickness, cex.lab = font_size, cex.axis = font_size)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  abline(h = 0, col = 'black', lty = 2, lwd = thickness)
  axis(2, col="black",lwd=2, cex.axis = axis_size)
  mtext(2,text = expression(v[Delta](x) - v(x)), line = 2, cex = font_size)
  axis(1, col="black", lwd=2, cex.axis = axis_size)
  mtext(1,text=paste("debt-to-income x",sep=""), line=2.5, cex = font_size)
  abline(v = x_bar, col = 'blue', lty = 4, lwd = thickness)
  abline(v = x_bar_nc, col = 'red', lty = 5, lwd = thickness)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # generate plot for debt price function
  y_min <- min(d,d_nc,d_c)
  y_max <- max(d,d_nc,d_c)
  pdf('poisson_d.pdf')
  par(mar=c(4, 4, 4, 4))
  plot(x = x, y = d, type='l', axes=F, ylim=c(y_min,y_max), 
       xlab="", ylab="", col="blue", main="", lwd = thickness, cex.lab = font_size, cex.axis = font_size)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  lines(x = x[-1], y = d_nc, col = 'red', lty = 2, lwd = thickness)
  lines(x = x[-1], y = d_c, col = 'purple', lty = 3, lwd = thickness)
  axis(2, col="black",lwd=2, cex.axis = axis_size)
  mtext(2,text = 'debt price', line = 2, cex = font_size)
  axis(1, col="black", lwd=2, cex.axis = axis_size)
  mtext(1,text=paste("debt-to-income x",sep=""), line=2.5, cex = font_size)
  abline(v = x_bar, col = 'blue', lty = 4, lwd = thickness)
  abline(v = x_bar_nc, col = 'red', lty = 5, lwd = thickness)
  abline(v = x_bar_c, col = 'purple', lty = 6, lwd = thickness)
  legend("bottomleft", legend = c('poisson','no-commitment','commitment'), 
         col = c('blue', 'red','purple'), bty="n", lty = c(1,2,3), cex = font_size, lwd = thickness)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # generate plot for expected issuance rate
  y_min <- min(exp_iss,g_nc)
  y_max <- 2
  pdf('poisson_g.pdf')
  par(mar=c(4, 4, 4, 4))
  plot(x = x, y = exp_iss, type='l', axes=F, ylim=c(y_min,y_max), 
       xlab="", ylab="", col="blue", main="", lwd = thickness, cex.lab = font_size, cex.axis = font_size)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  lines(x = x[-1], y = g_nc, col = 'red', lty = 2, lwd = thickness)
  lines(x = x, y = m*x, col = 'green', lty = 3, lwd = thickness)
  abline(v = x_a, col = 'lightgreen', lty = 3, lwd = thickness)
  axis(2, col="black",lwd=2, cex.axis = axis_size)
  mtext(2,text = 'expected issuance rate', line = 2, cex = font_size)
  axis(1, col="black", lwd=2, cex.axis = axis_size)
  mtext(1,text=paste("debt-to-income x",sep=""), line=2.5, cex = font_size)
  abline(v = x_bar, col = 'blue', lty = 4, lwd = thickness)
  abline(v = x_bar_nc, col = 'red', lty = 5, lwd = thickness)
  legend("topright", legend = c('poisson','no-commitment'), 
         col = c('blue', 'red'), bty="n", lty = c(1,2), cex = font_size, lwd = thickness)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # generate plot for 'jump-to'
  y_min <- min(jump_to,x)
  y_max <- max(jump_to,x)
  pdf('poisson_jump_to.pdf')
  par(mar=c(4, 4, 4, 4))
  plot(x = x, y = jump_to, type='l', axes=F, ylim=c(y_min,y_max), 
       xlab="", ylab="", col="blue", main="", lwd = thickness, cex.lab = font_size, cex.axis = font_size)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  lines(x = x, y = x, col = 'red', lty = 2, lwd = thickness)
  axis(2, col="black",lwd=2, cex.axis = axis_size)
  mtext(2,text = 'jump to', line = 2, cex = font_size)
  axis(1, col="black", lwd=2, cex.axis = axis_size)
  mtext(1,text=paste("debt-to-income x",sep=""), line=2.5, cex = font_size)
  abline(v = x_bar, col = 'blue', lty = 4, lwd = thickness)
  abline(v = x_bar_nc, col = 'red', lty = 5, lwd = thickness)
  legend("topright", legend = c('poisson','no-commitment'), 
         col = c('blue', 'red'), bty="n", lty = c(1,2), cex = font_size, lwd = thickness)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # generate plot for 'expected jump value'
  y_min <- min(jump_val,0)
  y_max <- max(jump_val,0)
  pdf('poisson_jump_val.pdf')
  par(mar=c(4, 4, 4, 4))
  plot(x = x, y = jump_val, type='l', axes=F, ylim=c(y_min,y_max), 
       xlab="", ylab="", col="blue", main="", lwd = thickness, cex.lab = font_size, cex.axis = font_size)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  abline(h = 0, col = 'red', lty = 2, lwd = thickness)
  axis(2, col="black",lwd=2, cex.axis = axis_size)
  mtext(2,text = 'expected jump value', line = 2, cex = font_size)
  axis(1, col="black", lwd=2, cex.axis = axis_size)
  mtext(1,text=paste("debt-to-income x",sep=""), line=2.5, cex = font_size)
  abline(v = x_bar, col = 'blue', lty = 4, lwd = thickness)
  abline(v = x_bar_nc, col = 'red', lty = 5, lwd = thickness)
  legend("topright", legend = c('poisson','no-commitment'), 
         col = c('blue', 'red'), bty="n", lty = c(1,2), cex = font_size, lwd = thickness)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # generate plot for error
  pdf('calc_errors.pdf')
  par(mar=c(4, 4, 4, 4))
  plot(x = 2:iter, y = err_t[-1], type='l', axes=F, xlab="", ylab="", col="blue", 
       main="", lwd = thickness, cex.lab = font_size, cex.axis = font_size, ylim = c(0,tol*100))
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  axis(2, col="black",lwd=2, cex.axis = axis_size)
  mtext(2,text = 'iteration error', line = 2, cex = font_size)
  axis(1, col="black", lwd=2, cex.axis = axis_size)
  mtext(1,text=paste("iteration number",sep=""), line=2.5, cex = font_size)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # compute stationary moments
  eq_object              <- list()
  eq_object[['x']]       <- x
  eq_object[['exp_iss']] <- exp_iss
  eq_object[['d']]       <- d
  eq_object[['v']]       <- v
  eq_object[['jump_to']] <- jump_to
  eq_object[['jump_val']]<- jump_val
  
  # smooth functions for output purposes
  eq_fn        <- list()
  for (func in names(eq_object)){
    eq_fn[[func]] <- splinefun(x, eq_object[[func]], method = "natural")
  }
  
  # other moments to report
  other_mom               <- c()
  other_mom['x_bar']      <- x_bar
  other_mom['v_0']        <- v_0
  other_mom['x_a']        <- x_a
  
  # store algo performance metrics
  perf_list               <- c()
  perf_list['nb_iter']    <- iter
  perf_list['error']      <- err
  
  # create output that will be returned by the program
  output        <- list(eq_fn, eq_object, other_mom, num_param, perf_list)
  output_names  <- c("eq_functions", "eq_objects", "other_moments", "num_param", "algo_perf")
  names(output) <- output_names
  
  # export data in data table
  out_df        <- as.data.frame(matrix(0, ncol = length(eq_object), nrow = n+1))
  names(out_df) <- names(eq_object)
  for (dummy in names(out_df)){
    out_df[,dummy] <- eq_object[[dummy]]
  }
  
  # export all policies and stationary distribution into csv file
  write.table(out_df, file = paste("policies.csv",sep=""), eol = "\n", na = "NA", fileEncoding = "",
              row.names = FALSE)
  return(output)
}

Poisson_commit_debt_root_fn <- function(param){
  # compute root eta for debt function with commitment and Poisson financing
  delta   <- param["delta"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  lambda  <- param["lambda"]
  r       <- param["r"]
  nu      <- param["nu"]
  
  # compute root of characteristic polynomial for debt HJB
  a       <- 0.5*sigma^2
  b       <- - (m + mu - nu * sigma - 0.5*sigma^2)
  c       <- - (r + m + lambda)
  discrim <- b^2 - 4*a*c
  eta     <- (- b + sqrt(discrim))/(2*a)
  
  return(eta)
}

Poisson_commit_debt_root_fn2 <- function(param){
  # compute roots eta_1 and eta_2 for debt function with commitment and Poisson financing
  delta   <- param["delta"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  lambda  <- param["lambda"]
  r       <- param["r"]
  nu      <- param["nu"]
  
  # compute root of characteristic polynomial for debt HJB
  a       <- 0.5*sigma^2
  b       <- - (m + mu - nu * sigma - 0.5*sigma^2)
  c       <- - (r + m + lambda)
  discrim <- b^2 - 4*a*c
  eta_1   <- (- b - sqrt(discrim))/(2*a)
  eta_2   <- (- b + sqrt(discrim))/(2*a)
  
  eta     <- c(eta_1,eta_2)
  names(eta) <- c("eta_1", "eta_2")
  return(eta)
}

Poisson_commit_value_root_fn <- function(param){
  # compute root xi for value function with commitment and Poisson financing
  delta   <- param["delta"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  lambda  <- param["lambda"]
  
  # compute root of characteristic polynomial for value function HJB
  a       <- 0.5*sigma^2
  b       <- - (m + mu + 0.5*sigma^2)
  c       <- - (delta + lambda - mu)
  discrim <- b^2 - 4*a*c
  xi      <- (- b + sqrt(discrim))/(2*a)
  
  return(xi)
}

Poisson_commit_value_root_fn2 <- function(param){
  # compute root xi_1 and xi_2 for value function with commitment and Poisson financing
  delta   <- param["delta"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  lambda  <- param["lambda"]
  
  # compute root of characteristic polynomial for value function HJB
  a       <- 0.5*sigma^2
  b       <- - (m + mu + 0.5*sigma^2)
  c       <- - (delta + lambda - mu)
  discrim <- b^2 - 4*a*c
  xi_2    <- (- b + sqrt(discrim))/(2*a)
  xi_1    <- (- b - sqrt(discrim))/(2*a)
  
  xi      <- c(xi_1,xi_2)
  names(xi) <- c("xi_1", "xi_2")
  return(xi)
}

#===================================================
# Value function Poisson opportunities to sell bonds
# commitment to always go back to debt-to-income x*
# default if jump to x* is worse than default
#===================================================

Poisson_commit_debt_px_fn2 <- function(param, x, x_star, x_bounds){
  # compute debt price for case where financing opportunities arrive at Poisson times 
  
  # load parameters
  delta     <- param["delta"]
  lambda    <- param["lambda"]
  kappa     <- param["kappa"]
  m         <- param["m"]
  mu        <- param["mu"]
  nu        <- param["nu"]
  sigma     <- param["sigma"]
  alpha     <- param["alpha"]
  theta     <- param["theta"]
  r         <- param["r"]
  
  # compute optimal boundaries with commitment
  x_bar     <- x_bounds['x_bar']
  x_bar_bar <- x_bounds['x_bar_bar']
  rho_star  <- x_star/x_bar
  rho_bar   <- x_bar/x_bar_bar
  
  # compute debt root
  debt_root <- Poisson_commit_debt_root_fn2(param)
  eta_1     <- debt_root['eta_1']
  eta_2     <- debt_root['eta_2']
  
  # set up system of 4 equations in 4 unknown to compute d*, d_bar, d1, d2
  A_mat     <- array(0, dim = c(4,4))
  B_vec     <- array(0, 4)
  
  # equation 1: equation linking d* and d_bar
  A_mat[1,1]<- 1
  A_mat[1,2]<- -(r+m+lambda)/(r+m+lambda*rho_star^eta_2)*rho_star^eta_2
  B_vec[1]  <- (kappa+m)*(1 - rho_star^eta_2)/(r+m+lambda*rho_star^eta_2)
  
  # equation 2: value matching at x = x_bar
  A_mat[2,2]<- -1
  A_mat[2,3]<- rho_bar^eta_1
  A_mat[2,4]<- rho_bar^eta_2
  B_vec[2]  <- -(kappa+m)/(r+m+lambda)
  
  # equation 3: smooth pasting at x = x_bar
  A_mat[3,1]<- eta_2*lambda/(r + m + lambda)
  A_mat[3,2]<- -eta_2
  A_mat[3,3]<- eta_1*rho_bar^eta_1
  A_mat[3,4]<- eta_2*rho_bar^eta_2
  B_vec[3]  <- -eta_2*(kappa+m)/(r + m + lambda)
  
  # equation 4: default boundary at x = x_bar_bar
  A_mat[4,3]<- 1
  A_mat[4,4]<- 1
  B_vec[4]  <- -(kappa+m)/(r + m + lambda)
  
  # compute vector of unknown d*, d_bar, d1, d2
  d_vec   <- solve(A_mat, B_vec)
  d_star  <- d_vec[1]
  d_bar   <- d_vec[2]
  d1      <- d_vec[3]
  d2      <- d_vec[4]
  
  # compute debt price for either x < x_bar or x > x_bar
  d_case_1  <- (kappa+m+lambda*d_star)/(r+m+lambda) - 
    ((kappa+m+lambda*d_star)/(r+m+lambda) - d_bar)*(x/x_bar)^eta_2
  d_case_2  <- (kappa+m)/(r+m+lambda) + 
    d1*(x/x_bar_bar)^eta_1 + d2*(x/x_bar_bar)^eta_2
  
  # compute d(x)
  d     <- d_case_1*(x <= x_bar) + d_case_2*(x > x_bar)*(x <= x_bar_bar) + 0*(x > x_bar_bar)
  return(d)
}

Poisson_commit_value_fn2 <- function(param, x, x_star, x_bounds){
  # compute value function for case where financing opportunities arrive at Poisson times 
  
  # load parameters
  delta     <- param["delta"]
  lambda    <- param["lambda"]
  kappa     <- param["kappa"]
  m         <- param["m"]
  mu        <- param["mu"]
  nu        <- param["nu"]
  sigma     <- param["sigma"]
  alpha     <- param["alpha"]
  theta     <- param["theta"]
  r         <- param["r"]
  
  # compute optimal boundaries with commitment
  x_bar     <- x_bounds['x_bar']
  x_bar_bar <- x_bounds['x_bar_bar']
  rho_star  <- x_star/x_bar
  rho_bar   <- x_bar/x_bar_bar
  
  # compute debt root
  debt_root <- Poisson_commit_debt_root_fn2(param)
  eta_1     <- debt_root['eta_1']
  eta_2     <- debt_root['eta_2']
  
  # compute value root
  value_root<- Poisson_commit_value_root_fn2(param)
  xi_1      <- value_root['xi_1']
  xi_2      <- value_root['xi_2']
  
  # first, compute debt constants d*, d_bar, d1, d2
  # set up system of 4 equations in 4 unknown to compute d*, d_bar, d1, d2
  A_mat     <- array(0, dim = c(4,4))
  B_vec     <- array(0, 4)
  
  # equation 1: equation linking d* and d_bar
  A_mat[1,1]<- 1
  A_mat[1,2]<- -(r+m+lambda)/(r+m+lambda*rho_star^eta_2)*rho_star^eta_2
  B_vec[1]  <- (kappa+m)*(1 - rho_star^eta_2)/(r+m+lambda*rho_star^eta_2)
  
  # equation 2: value matching at x = x_bar
  A_mat[2,2]<- -1
  A_mat[2,3]<- rho_bar^eta_1
  A_mat[2,4]<- rho_bar^eta_2
  B_vec[2]  <- -(kappa+m)/(r+m+lambda)
  
  # equation 3: smooth pasting at x = x_bar
  A_mat[3,1]<- eta_2*lambda/(r + m + lambda)
  A_mat[3,2]<- -eta_2
  A_mat[3,3]<- eta_1*rho_bar^eta_1
  A_mat[3,4]<- eta_2*rho_bar^eta_2
  B_vec[3]  <- -eta_2*(kappa+m)/(r + m + lambda)
  
  # equation 4: default boundary at x = x_bar_bar
  A_mat[4,3]<- 1
  A_mat[4,4]<- 1
  B_vec[4]  <- -(kappa+m)/(r + m + lambda)
  
  # compute vector of unknown d*, d_bar, d1, d2
  d_vec   <- solve(A_mat, B_vec)
  d_star  <- d_vec[1]
  d_bar   <- d_vec[2]
  d1      <- d_vec[3]
  d2      <- d_vec[4]
  
  # then, compute value function constants v*, v_bar, v1, v2
  # set up system of 4 equations in 4 unknown to compute v*, v_bar, v1, v2
  C_mat     <- array(0, dim = c(4,4))
  D_vec     <- array(0, 4)
  
  # equation 1: equation linking v* and v_bar
  C_mat[1,1]<- 1 - lambda*(1 - rho_star^xi_2)/(delta + lambda - mu)
  C_mat[1,2]<- -rho_star^xi_2
  D_vec[1]  <- (1 + lambda*x_star*d_star)*(1 - rho_star^xi_2)/(delta + lambda - mu) - 
    x_star*(kappa+m+lambda*d_star)/(delta+m+lambda)*(1-rho_star^(xi_2-1))
  
  # equation 2: value matching at x = x_bar
  C_mat[2,2]<- -1
  C_mat[2,3]<- rho_bar^xi_1
  C_mat[2,4]<- rho_bar^xi_2
  D_vec[2]  <- (kappa+m)/(delta+m+lambda)*x_bar - 1/(delta+lambda-mu)
  
  # equation 3: smooth pasting at x = x_bar
  C_mat[3,1]<- xi_2*lambda/(delta + lambda - mu)
  C_mat[3,2]<- -xi_2
  C_mat[3,3]<- xi_1*rho_bar^xi_1
  C_mat[3,4]<- xi_2*rho_bar^xi_2
  D_vec[3]  <- (kappa+m)/(delta + m + lambda)*x_bar - 
    xi_2*(1 + lambda*x_star*d_star)/(delta + lambda - mu) + 
    (xi_2-1)*(kappa+m+lambda*d_star)/(delta + m + lambda)*x_bar
  
  # equation 4: default boundary at x = x_bar_bar
  C_mat[4,3]<- 1
  C_mat[4,4]<- 1
  D_vec[4]  <- (kappa+m)/(delta + m + lambda)*x_bar_bar - 1/(delta+lambda-mu)
  
  # compute vector of unknown v*, v_bar, v1, v2
  v_vec     <- solve(C_mat, D_vec)
  v_star    <- v_vec[1]
  v_bar     <- v_vec[2]
  v1        <- v_vec[3]
  v2        <- v_vec[4]
  
  # compute value function for either x < x_bar or x > x_bar
  v_case_1  <- (1+lambda*(v_star+d_star*x_star))/(delta+lambda-mu)*(1-(x/x_bar)^xi_2) - 
    x*(kappa+m+lambda*d_star)/(delta+m+lambda)*(1-(x/x_bar)^(xi_2-1)) + 
    v_bar*(x/x_bar)^xi_2
  v_case_2  <- 1/(delta+lambda-mu) - x*(kappa+m)/(delta+m+lambda) + 
    v1*(x/x_bar_bar)^xi_1 + v2*(x/x_bar_bar)^xi_2 
  
  # compute v(x)
  v     <- v_case_1*(x <= x_bar) + v_case_2*(x > x_bar)*(x <= x_bar_bar) + 0*(x > x_bar_bar)
  return(v)
}

Poisson_commit_cutoff_fn2 <- function(param, x_star){
  # compute optimal boundaries for case where financing opportunities arrive at Poisson times
  # x_star is given
  # compute x_bar and x_bar_bar
  
  # numerical parameters
  eps       <- 1E-4     # threshold for whether solution exists or not
  nb        <- 20       # number of points at which to compute implicit equations
  min_dist  <- .25      # minimum distance between x*, x_bar and x_bar_bar
  min_r_cond<- 1E-200   # minimum conditioning number
  
  # load parameters
  delta     <- param["delta"]
  lambda    <- param["lambda"]
  kappa     <- param["kappa"]
  m         <- param["m"]
  mu        <- param["mu"]
  nu        <- param["nu"]
  sigma     <- param["sigma"]
  alpha     <- param["alpha"]
  theta     <- param["theta"]
  r         <- param["r"]
  
  # compute debt root
  debt_root <- Poisson_commit_debt_root_fn2(param)
  eta_1     <- debt_root['eta_1']
  eta_2     <- debt_root['eta_2']
  
  # compute value root
  value_root<- Poisson_commit_value_root_fn2(param)
  xi_1      <- value_root['xi_1']
  xi_2      <- value_root['xi_2']
  
  # given x_bar and x_bar_bar, compute the 2 non-linear equations that need to be solved
  implicit_fn <- function(x_bar, x_bar_bar){
    # given x_bar and x_bar_bar, compute constants for debt price and value function
    rho_star  <- x_star/x_bar
    rho_bar   <- x_bar/x_bar_bar
    
    # first, compute debt constants d*, d_bar, d1, d2
    # set up system of 4 equations in 4 unknown to compute d*, d_bar, d1, d2
    A_mat     <- array(0, dim = c(4,4))
    B_vec     <- array(0, 4)
    
    # equation 1: equation linking d* and d_bar
    A_mat[1,1]<- 1
    A_mat[1,2]<- -(r+m+lambda)/(r+m+lambda*rho_star^eta_2)*rho_star^eta_2
    B_vec[1]  <- (kappa+m)*(1 - rho_star^eta_2)/(r+m+lambda*rho_star^eta_2)
    
    # equation 2: value matching at x = x_bar
    A_mat[2,2]<- -1
    A_mat[2,3]<- rho_bar^eta_1
    A_mat[2,4]<- rho_bar^eta_2
    B_vec[2]  <- -(kappa+m)/(r+m+lambda)
    
    # equation 3: smooth pasting at x = x_bar
    A_mat[3,1]<- eta_2*lambda/(r + m + lambda)
    A_mat[3,2]<- -eta_2
    A_mat[3,3]<- eta_1*rho_bar^eta_1
    A_mat[3,4]<- eta_2*rho_bar^eta_2
    B_vec[3]  <- -eta_2*(kappa+m)/(r + m + lambda)
    
    # equation 4: default boundary at x = x_bar_bar
    A_mat[4,3]<- 1
    A_mat[4,4]<- 1
    B_vec[4]  <- -(kappa+m)/(r + m + lambda)
    
    # compute vector of unknown d*, d_bar, d1, d2
    d_vec   <- solve(A_mat, B_vec, tol = min_r_cond)
    d_star  <- d_vec[1]
    d_bar   <- d_vec[2]
    d1      <- d_vec[3]
    d2      <- d_vec[4]
    
    # then, compute value function constants v*, v_bar, v1, v2
    # set up system of 4 equations in 4 unknown to compute v*, v_bar, v1, v2
    C_mat     <- array(0, dim = c(4,4))
    D_vec     <- array(0, 4)
    
    # equation 1: equation linking v* and v_bar
    C_mat[1,1]<- 1 - lambda*(1 - rho_star^xi_2)/(delta + lambda - mu)
    C_mat[1,2]<- -rho_star^xi_2
    D_vec[1]  <- (1 + lambda*x_star*d_star)*(1 - rho_star^xi_2)/(delta + lambda - mu) - 
      x_star*(kappa+m+lambda*d_star)/(delta+m+lambda)*(1-rho_star^(xi_2-1))
    
    # equation 2: value matching at x = x_bar
    C_mat[2,2]<- -1
    C_mat[2,3]<- rho_bar^xi_1
    C_mat[2,4]<- rho_bar^xi_2
    D_vec[2]  <- (kappa+m)/(delta+m+lambda)*x_bar - 1/(delta+lambda-mu)
    
    # equation 3: smooth pasting at x = x_bar
    C_mat[3,1]<- xi_2*lambda/(delta + lambda - mu)
    C_mat[3,2]<- -xi_2
    C_mat[3,3]<- xi_1*rho_bar^xi_1
    C_mat[3,4]<- xi_2*rho_bar^xi_2
    D_vec[3]  <- (kappa+m)/(delta + m + lambda)*x_bar - 
      xi_2*(1 + lambda*x_star*d_star)/(delta + lambda - mu) + 
      (xi_2-1)*(kappa+m+lambda*d_star)/(delta + m + lambda)*x_bar
    
    # equation 4: default boundary at x = x_bar_bar
    C_mat[4,3]<- 1
    C_mat[4,4]<- 1
    D_vec[4]  <- (kappa+m)/(delta + m + lambda)*x_bar_bar - 1/(delta+lambda-mu)
    
    # compute vector of unknown v*, v_bar, v1, v2
    v_vec     <- solve(C_mat, D_vec, tol = min_r_cond)
    v_star    <- v_vec[1]
    v_bar     <- v_vec[2]
    v1        <- v_vec[3]
    v2        <- v_vec[4]
    
    # compute implicit function 0 = v_star + (x_star - x_bar)*d_star
    condition_1 <- v_star + (x_star - x_bar)*d_star
    
    # compute implicit function v'(x_bar_bar) = 0
    condition_2 <- - (kappa+m)/(delta+m+lambda)*x_bar_bar + v1*xi_1 + v2*xi_2 
    
    # conditions
    conditions  <- c(condition_1,condition_2)
    return(conditions)
  }

  # given x_bar and x_bar_bar, compute the 2 non-linear equations that need to be solved
  # return only the first of these 2 conditions
  implicit_fn1 <- function(x_bar, x_bar_bar){
    # given x_bar and x_bar_bar, compute constants for debt price and value function
    return(implicit_fn(x_bar, x_bar_bar)[1])
  }
  
  # given x_bar and x_bar_bar, compute the 2 non-linear equations that need to be solved
  # return only the second of these 2 conditions
  implicit_fn2 <- function(x_bar, x_bar_bar){
    # given x_bar and x_bar_bar, compute constants for debt price and value function
    return(implicit_fn(x_bar, x_bar_bar)[2])
  }
  
  # initial search for x_bar_bar -- the default boundary
  x_bb_test     <- seq(from = x_star + min_dist, to = 1.5*no_commit_cutoff_fn(param), length.out = nb)
  out_x_bb_test <- array(NA,nb)
  step_size     <- mean(diff(x_bb_test))
  iter          <- 0
  
  # loop
  while (step_size > .0025){
    for (j in (1:nb)){
      # set x_bar_bar = x_bb, and search for x_bar
      x_bb         <- x_bb_test[j]
      
      # look for x_b given x_bb, using a grid initially
      x_b_test     <- seq(from = x_star+min_dist, to = x_bb-min_dist, length.out = nb)
      out_x_b_test <- apply(X = as.matrix(1:nb), 1, 
                            FUN = function(i) implicit_fn1(x_b_test[i], x_bar_bar = x_bb))
      
      # check that min is negative and max is positive
      min_test_x_b <- min(out_x_b_test)
      max_test_x_b <- max(out_x_b_test)
      sign_end_x_b <- sign(out_x_b_test[nb])
      
      # check if there is at least one solution
      if ((min_test_x_b>-eps)|(max_test_x_b<eps)){
        # case where there is no solution
        x_b    <- NA
      } else {
        # case where there is at least one solution -- in this case, pick highest solution
        x_low  <- x_b_test[max(which(sign(out_x_b_test)!=sign_end_x_b))]
        x_high <- x_bb - min_dist
        x_b    <- uniroot(f = implicit_fn1, x_bar_bar = x_bb, interval = c(x_low,x_high))$root
      }
      
      if (!is.na(x_b)){
        out_x_bb_test[j] <- implicit_fn2(x_b,x_bb)
      } else {
        out_x_bb_test[j] <- NA
      }
    }
    
    # look at min and max for out_x_bb_test, removing all NAs
    x_bb_test     <- x_bb_test[!is.na(out_x_bb_test)]
    out_x_bb_test <- out_x_bb_test[!is.na(out_x_bb_test)]
    
    # if length of out_x_bb_test is below 1, no solution
    if (length(x_bb_test) < 2){
      # no solution
      x_bb <- NA
      break
    } else {
      # check min and max
      min_test_x_bb <- min(out_x_bb_test)
      max_test_x_bb <- max(out_x_bb_test)
      sign_end_x_bb <- sign(out_x_bb_test[length(out_x_bb_test)])
      
      # if min and max have same sign, then x_bb = NA
      if ((min_test_x_bb>-eps)|(max_test_x_bb<eps)){
        # case where there is no solution
        x_bb <- NA
        break
      } else {
        # case where there is a solution
        i_interim <- max(which(sign(out_x_bb_test)!=sign_end_x_bb))
        x_bb      <- .5*(x_bb_test[i_interim] + x_bb_test[i_interim+1])
        # update vector of x_bb to test
        x_bb_test     <- seq(from = x_bb - 2*step_size, to = x_bb + 2*step_size, length.out = nb)
        out_x_bb_test <- array(NA,nb)
        step_size     <- mean(diff(x_bb_test))
      }
    }
    
    # next iteration -- update counter
    iter          <- iter+1
  }
  
  # for x_bb found, update x_b
  if (!is.na(x_bb)){
    x_b_test     <- seq(from = x_star + min_dist, to = x_bb - min_dist, length.out = nb)
    out_x_b_test <- apply(X = as.matrix(1:nb), 1, 
                          FUN = function(i) implicit_fn1(x_b_test[i],x_bar_bar = x_bb))
    
    # check that min is negative and max is positive
    min_test_x_b <- min(out_x_b_test)
    max_test_x_b <- max(out_x_b_test)
    sign_end_x_b <- sign(out_x_b_test[nb])
    
    # check if there is at least one solution
    if ((min_test_x_b>-eps)|(max_test_x_b<eps)){
      # case where there is no solution
      x_b    <- NA
    } else {
      # case where there is at least one solution -- in this case, pick highest solution
      x_low  <- x_b_test[max(which(sign(out_x_b_test)!=sign_end_x_b))]
      x_high <- x_bb - min_dist
      x_b    <- uniroot(f = implicit_fn1, x_bar_bar = x_bb, interval = c(x_low,x_high))$root
    }    
    
    # update objectives
    objectives        <- implicit_fn(x_b, x_bb)
  } else {
    x_b        <- NA
    x_bb       <- NA
    objectives <- c(NA,NA)
  }
  
  # return boundaries
  boundaries        <- c(x_b, x_bb)
  names(boundaries) <- c("x_bar", "x_bar_bar")
  
  return(boundaries)
}

Poisson_commit_optim_leverage_fn2 <- function(param){
  # function that finds optimal x* for model with debt-to-income commitment
  
  # grid search parameters
  nb         <- 100
  eps        <- 1E-3
  low_bound  <- no_commit_cutoff_fn(param)/10
  high_bound <- no_commit_cutoff_fn(param)
  dist_bound <- high_bound - low_bound
  
  dummy      <- 1
  while (dist_bound > eps){
    # print iteration
    print(paste0('iteration: ',dummy))
    
    # compute initial x_star to test
    test       <- seq(from = low_bound, to = high_bound, length.out = nb)
    objective  <- rep(NA,nb)
    cutoff     <- rep(NA,nb)
    
    # compute objective for each possible x*
    for (i in (1:nb)){
      prout     <- Poisson_commit_cutoff_fn2(param, x_star = test[i])
      cutoff[i] <- prout['x_bar_bar']
      if (!is.na(cutoff[i])) {
        objective[i] <- Poisson_commit_value_fn2(param, x = test[i], x_star = test[i], x_bounds = prout) + 
          test[i] * Poisson_commit_debt_px_fn2(param, x = test[i], x_star = test[i], x_bounds = prout)
      }
    }
    
    # select only those tests that do not return NA
    test      <- test[!is.na(cutoff)]
    objective <- objective[!is.na(cutoff)]
    
    # find best x*
    x_opt <- test[which.max(objective)]
    
    # update bounds of search
    low_bound  <- x_opt - dist_bound/20
    high_bound <- x_opt + dist_bound/20
    dist_bound <- dist_bound/10
    dummy      <- dummy + 1
  }
  
  # return optimal initial debt-to-income
  return(x_opt)
}

plot_Poisson_commit_fn2 <- function(param, x_star, plotting_param){
  # function that plots all useful plots
  # plot functions for Poisson arrival of debt issuance opportunities
  # commitment to a particular capital structure
  
  # plotting parameters
  n              <- plotting_param['n']
  thickness      <- plotting_param['thick']
  font_size      <- plotting_param['font_size']
  axis_size      <- plotting_param['axis_size']
  
  # compute boundaries
  x_bounds <- Poisson_commit_cutoff_fn2(param, x_star)
  x_bar    <- x_bounds['x_bar']
  x_bar_bar<- x_bounds['x_bar_bar']
  
  # create grid
  x_min    <- .01
  x_max    <- x_bar_bar
  x        <- seq(from = x_min, to = x_max, length.out = n)
  
  # compute v and d at x* 
  v_star   <- Poisson_commit_value_fn2(param, x_star, x_star, x_bounds)
  d_star   <- Poisson_commit_debt_px_fn2(param, x_star, x_star, x_bounds)
  
  # construct grid of output vectors
  out_list                <- list()
  out_list[['debt']]      <- Poisson_commit_debt_px_fn2(param, x, x_star, x_bounds)
  out_list[['value']]     <- Poisson_commit_value_fn2(param, x, x_star, x_bounds)
  out_list[['check']]     <- v_star+(x_star-x)*d_star
  
  # name plots
  legend <- c('debt price', 'value function', 'check')
  names(legend) <- names(out_list)
  
  # plot all functions of interest
  for (plot in names(out_list)){
    y_min <- min(out_list[[plot]], na.rm = TRUE)
    y_max <- max(out_list[[plot]], na.rm = TRUE)
    pdf(paste0(gsub(" ", "_", legend[plot]),"_Poisson_model_commit.pdf"))
    par(mar=c(4, 4, 4, 4))
    plot(x = x, y = out_list[[plot]], type='l', axes=F, xlim=c(x_min,x_max), ylim=c(y_min,y_max), lty = 1, 
         xlab="", ylab="",col="blue", main="", lwd = thickness, cex.lab = font_size, cex.axis = font_size)
    grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
    u <- par('usr')
    if (plot == 'check') abline(h = 0, col = 'green', lty = 3, lwd = thickness)
    if (plot == 'value') abline(h = 0, col = 'black', lty = 3, lwd = thickness)
    if (plot == 'value') lines(x = x, y = out_list[['check']], col = 'lightblue', lty = 2, lwd = thickness)
    abline(v = x_bar, col = 'lightblue', lty = 4, lwd = thickness)
    abline(v = x_bar_bar, col = 'red', lty = 5, lwd = thickness)
    abline(v = x_star, col = 'purple', lty = 6, lwd = thickness)
    axis(2, col="black",lwd=2, cex.axis = axis_size)
    mtext(2,text = legend[plot], line = 2, cex = font_size)
    axis(1, col="black", lwd=2, cex.axis = axis_size)
    mtext(1,text=paste("debt-to-income x",sep=""),line=2.5, cex = font_size)
    if (plot == 'value') legend("topright", legend = c('value function', 'value upon jump to x*'), col = c('blue', 'lightblue'), bty="n", cex = font_size, lwd = thickness, lty = 1:2)
    text(x = x_star - .04*(u[2]-u[1]), y = u[3] + .6*(u[4]-u[3]), "x*", cex = font_size, col = 'black') 
    text(x = x_bar - .04*(u[2]-u[1]), y = u[3] + .6*(u[4]-u[3]), expression(bar(x)), cex = font_size, col = 'black') 
    text(x = x_bar_bar - .04*(u[2]-u[1]), y = u[3] + .6*(u[4]-u[3]), expression(bar(bar(x))), cex = font_size, col = 'black') 
    box(which = "plot", bty = "l", lwd = thickness)
    dev.off()
  }
}

#=======================================
# Markov Switching Equilibrium Functions
#=======================================

markov_switching_root_fn <- function(param){
  # compute root eta for debt price with markov switching model
  # load parameters
  r       <- param["r"]
  lambda_c<- param["lambda_c"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  
  # compute root of characteristic polynomial
  a       <- 0.5*sigma^2
  b       <- - (m + mu - 0.5*sigma^2)
  c       <- - (r + m + lambda_c)
  discrim <- b^2 - 4*a*c
  eta     <- (- b + sqrt(discrim))/(2*a)
  
  return(eta)
}

markov_switching_debt_px_u_fn <- function(param, x){
  # compute debt price in state "u" for markov switching equilibrium
  # no ability to issue debt in state "c"
  # full ability to issue debt in state "u"
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  lambda_c<- param["lambda_c"]
  
  # compute roots of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  eta     <- markov_switching_root_fn(param)
  x_bar   <- no_commit_cutoff_fn(param)
  
  # compute debt price in state "u"
  d_rf    <- (kappa + m)/(delta + m)
  debt_u  <- d_rf * (1 -(x/x_bar)^(xi-1))
  
  return(debt_u)
}

markov_switching_debt_px_c_fn <- function(param, x){
  # compute debt price in state "c" for markov switching equilibrium
  # no ability to issue debt in state "c"
  # full ability to issue debt in state "u"
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  lambda_c<- param["lambda_c"]
  
  # compute roots of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  eta     <- markov_switching_root_fn(param)
  x_bar   <- no_commit_cutoff_fn(param)
  
  # compute debt price in state "c"
  d_rf    <- (kappa + m)/(delta + m)
  debt_c  <- (kappa + m)/(r + m + lambda_c) * 
    (1 + lambda_c/(delta+m)) + 
    lambda_c/(delta - (r+lambda_c))*d_rf*(x/x_bar)^(xi-1) - 
    ((kappa + m)/(r + m + lambda_c) * (1 + lambda_c/(delta+m)) + 
       lambda_c/(delta - (r+lambda_c))*d_rf)*(x/x_bar)^eta
  
  return(debt_c)
}

markov_switching_g_fn <- function(param, x){
  # compute issuance policy for markov switching equilibrium
  # no ability to issue debt in state "c"
  # full ability to issue debt in state "u"
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  lambda_c<- param["lambda_c"]
  lambda_u<- param["lambda_u"]
  
  # compute roots of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  eta     <- markov_switching_root_fn(param)
  x_bar   <- no_commit_cutoff_fn(param)
  
  # compute issuance strategy in unconstrained state
  d_rf    <- (kappa + m)/(delta + m)
  d_c     <- markov_switching_debt_px_c_fn(param, x)
  d_u     <- d_rf*(1 - (x/x_bar)^(xi-1))
  d_u_p   <- d_rf*(xi-1)/x_bar*(x/x_bar)^(xi-2)
  
  g    <- ((delta - (r+lambda_u))*d_u + lambda_u*d_c)/d_u_p
  
  return(g)
}

markov_switching_hat_value_fn <- function(param, num_param){
  # compute value function for citizens in case of Markov switching equilibrium
  
  # loan numerical parameters
  n_x       <- num_param['n_x']
  dt        <- num_param['dt']
  max_iter  <- num_param['max_iter']
  tol       <- num_param['tol']
  
  # load parameters
  delta     <- param["delta"]
  hat_delta <- param["hat_delta"]
  kappa     <- param["kappa"]
  m         <- param["m"]
  mu        <- param["mu"]
  sigma     <- param["sigma"]
  alpha     <- param["alpha"]
  theta     <- param["theta"]
  r         <- param["r"]
  lambda_c  <- param["lambda_c"]
  lambda_u  <- param["lambda_u"]
  
  # compute default boundary without commitment
  x_bar    <- no_commit_cutoff_fn(param)
  
  # set x_min so that we do not need to deal with singularity
  x_min    <- .01
  
  # create grid
  x        <- seq(from = x_min, to = x_bar, length.out = n_x)
  dx       <- x[2]-x[1]
  
  # compute debt price, consumption and issuance policy, as well as initial guess for hat v
  g_u     <- markov_switching_g_fn(param, x)
  g_c     <- rep(0,n_x)
  d_u     <- markov_switching_debt_px_u_fn(param, x)
  d_c     <- markov_switching_debt_px_c_fn(param, x)
  cons_u  <- 1 - (kappa+m)*x + g_u*d_u
  cons_c  <- 1 - (kappa+m)*x + g_c*d_c
  
  # guess initial value in both cases
  orig_v  <- no_commit_value_fn(param, x)
  v_n     <- c(orig_v,orig_v)
  
  # compute drift rate (invariant across iterations)
  drift_u <- g_u - (mu + m)*x
  drift_c <- g_c - (mu + m)*x
  
  # construct matrix A
  A        <- Matrix(0, nrow = 2*n_x, ncol = 2*n_x, sparse = TRUE)
  
  # diagonal elements for A
  locs     <- cbind((1:n_x),(1:n_x))
  A[locs]  <- -(lambda_u + abs(drift_u)/dx + (sigma*x/dx)^2)
  locs     <- cbind(((n_x+1):(2*n_x)),((n_x+1):(2*n_x)))
  A[locs]  <- -(lambda_c + abs(drift_c)/dx + (sigma*x/dx)^2)
  
  # lower diagonal elements -- probability to go to lower x
  locs     <- cbind((2:n_x),(1:(n_x-1)))
  A[locs]  <- -pmin(0,drift_u[-1])/dx + .5*(sigma*x[-1]/dx)^2
  locs     <- cbind(((n_x+2):(2*n_x)),((n_x+1):(2*n_x-1)))
  A[locs]  <- -pmin(0,drift_c[-1])/dx + .5*(sigma*x[-1]/dx)^2
  
  # upper diagonal elements -- probability to go to higher x
  locs     <- cbind((1:(n_x-1)),(2:n_x))
  A[locs]  <- pmax(0,drift_u[-n_x])/dx + .5*(sigma*x[-n_x]/dx)^2
  locs     <- cbind(((n_x+1):(2*n_x-1)),((n_x+2):(2*n_x)))
  A[locs]  <- pmax(0,drift_c[-n_x])/dx + .5*(sigma*x[-n_x]/dx)^2
  
  # transition into different markov state
  locs     <- cbind((1:n_x),((n_x+1):(2*n_x)))
  A[locs]  <- lambda_u
  locs     <- cbind(((n_x+1):(2*n_x)),(1:n_x))
  A[locs]  <- lambda_c
  
  # assume reflection at x = 0
  A[1,1]             <- -sum(A[1,-1])
  A[(n_x+1),(n_x+1)] <- -sum(A[(n_x+1),-(n_x+1)])
  A[n_x,n_x]         <- -sum(A[n_x,-n_x])
  A[(2*n_x),(2*n_x)] <- -sum(A[(2*n_x),-(2*n_x)])
  
  # construct B matrix and adjust for what happens at x = x_bar
  B                  <- (1/dt + hat_delta - mu)*Diagonal((2*n_x)) - A
  idx_reentry        <- which.min(abs(x-theta*x_bar))
  B[n_x,]            <- 0
  B[n_x,n_x]         <- 1
  B[n_x,idx_reentry] <- -alpha
  B[(2*n_x),]                  <- 0
  B[(2*n_x),(2*n_x)]           <- 1
  B[(2*n_x),(n_x+idx_reentry)] <- -alpha
  
  # constant part of the right hand side + take care of boundary conditions
  RHS_cst <- c(cons_u,cons_c)
  
  # initialization
  err     <- 1      # initialize error for value function
  err_vec <- c(1)
  iter    <- 1      # initialize iteration #
  # compute citizens' welfare value
  while ((err>tol)&(iter<max_iter)){
    # if error greater than tolerance
    v   <- v_n
    
    # enter right hand side, including correction for boundary conditions
    RHS        <- RHS_cst + v/dt
    RHS[n_x]   <- 0
    RHS[2*n_x] <- 0
    
    # solve linear system
    v_n  <- c(array(solve(B,RHS)))
    
    # compute error
    err      <- max(abs(v_n-v)/dt)
    err_vec  <- c(err_vec,err)
    
    # update iteration #
    iter     <- iter + 1
  }
  
  # interpolate at zero
  hat_v_0_u  <- v_n[1] - x_min*(v_n[2]-v_n[1])/dx
  hat_v_0_c  <- v_n[(n_x+1)] - x_min*(v_n[(n_x+2)]-v_n[(n_x+1)])/dx
  x          <- c(0,x)
  v_u        <- c(hat_v_0_u,v_n[1:n_x])
  v_c        <- c(hat_v_0_c,v_n[(n_x+1):(2*n_x)])
  
  # smooth function
  hat_v_u_fn <- approxfun(x, v_u, method = "linear")
  hat_v_c_fn <- approxfun(x, v_c, method = "linear")
  
  # output
  output        <- list(hat_v_u_fn,hat_v_0_u,hat_v_c_fn,hat_v_0_c,iter)
  names(output) <- c("hat_v_u", "hat_v0_u", "hat_v_c", "hat_v0_c", "n_iter")
  
  return(output)
}

#===================================
# Debt Ceiling Equilibrium Functions
#===================================

issuance_constraint_ratio_limit_fn <- function(param){
  # compute the minimum x*/bar x that is needed for a smooth equilibrium to exist
  # in which there is no issuance if debt-to-income is > bar x
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  r       <- param["r"]
  
  # compute roots of characteristic polynomial
  xi      <- no_commit_root_fn(param)
  eta_1   <- commit_debt_root2_fn(param)['eta_1']
  eta_2   <- commit_debt_root2_fn(param)['eta_2']
  
  # compute risk-free debt prices
  d_rf_delta <- (kappa+m)/(delta+m)
  d_rf_r     <- (kappa+m)/(r+m)
  
  # create implicit function whose root we need to find
  rho_fn <- function(rho){
    prout <- (eta_2*rho^(-eta_1) - eta_1*rho^(-eta_2))*
      (d_rf_r - d_rf_delta*(1 - rho^(xi-1))) - d_rf_r*(eta_2-eta_1)
    return(prout)
  }
  
  # find root of implicit function
  rho_lim <- uniroot(rho_fn, lower = .01, upper = .99, maxiter = 1000)$root
  return(rho_lim)
}

intermed_x_star_constraint_fn <- function(param){
  # compute the minimum x* above which the intermediate region exists
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  r       <- param["r"]
  mu      <- param["mu"]
  
  # step of grid search
  step    <- .001
  
  # first, find the x* below which the smooth equilibrium ceases to exist
  x_star_min <- issuance_constraint_ratio_limit_fn(param)*no_commit_cutoff_fn(param)
  
  # compute roots of characteristic polynomials
  xi_1    <- no_commit_root2_fn(param)['xi_1']
  xi_2    <- no_commit_root2_fn(param)['xi_2']
  eta_1   <- commit_debt_root2_fn(param)['eta_1']
  eta_2   <- commit_debt_root2_fn(param)['eta_2']
  
  # initialize search
  x_star  <- x_star_min - step
  flag    <- 1
  
  while(flag == 1){
    # update x_star
    x_star_n <- x_star
    
    # compute default boundary in reflecting equilibrium
    x_bar   <- constrained_reflecting_cutoff_fn(param,x_star_n)
    rho     <- x_star_n/x_bar
    
    # compute risk-free debt prices
    d_rf_delta <- (kappa+m)/(delta+m)
    d_rf_r     <- (kappa+m)/(r+m)
    
    # constants of integration in constrained region
    v_1     <- x_bar * (d_rf_delta*(1-xi_2*rho^(xi_2-1)) + 
                          xi_2*rho^(xi_2-1) / (delta-mu)/x_bar - 
                          d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
      (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
    v_2     <- x_bar * (d_rf_delta*(xi_1*rho^(xi_1-1) - 1) - 
                          xi_1*rho^(xi_1-1) / (delta-mu)/x_bar + 
                          d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
      (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
    
    # compute potential x_hat
    x_hat         <- xi_2*x_star_n/(1-xi_2)*
      ((1-xi_1)*v_1*rho^xi_1 + (1-xi_2)*v_2*rho^xi_2)/
      (xi_1*v_1*rho^xi_1 + xi_2*v_2*rho^xi_2)
    
    # update flag
    flag <- (x_hat>0)*(x_hat<x_star_n)
    
    # update x_star
    x_star <- x_star_n - step
  }
  
  return(x_star_n)
}

constrained_smooth_debt_px_fn <- function(param, x_star, x){
  # compute debt price for case with constrained for x > x*
  # assume smooth equilibrium
  delta   <- param["delta"]
  r       <- param["r"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  
  # compute roots of characteristic polynomials
  xi      <- no_commit_root_fn(param)
  eta_1   <- commit_debt_root2_fn(param)['eta_1']
  eta_2   <- commit_debt_root2_fn(param)['eta_2']
  x_bar   <- no_commit_cutoff_fn(param)
  rho     <- x_star/x_bar
  
  # compute risk-free debt prices
  d_rf_delta <- (kappa+m)/(delta+m)
  d_rf_r     <- (kappa+m)/(r+m)
  
  # compute debt price when unconstrained
  debt_u  <- d_rf_delta * (1 - (x/x_bar)^(xi - 1))
  
  # constants of integration in constrained region
  d_1     <- (d_rf_delta*(1 - rho^(xi-1)) - d_rf_r*(1-rho^eta_2))/(rho^eta_1 - rho^eta_2)
  d_2     <- (-d_rf_delta*(1 - rho^(xi-1)) + d_rf_r*(1-rho^eta_1))/(rho^eta_1 - rho^eta_2)
  
  # compute debt price when constrained
  debt_c  <- d_rf_r + d_1*(x/x_bar)^eta_1 + d_2*(x/x_bar)^eta_2
  
  # compute debt price
  debt    <- (x < x_star)*debt_u + (x >= x_star)*debt_c
  
  return(debt)
}

constrained_reflecting_debt_px_fn <- function(param, x_star, x){
  # compute debt price for case with constrained for x > x*
  # assume reflecting equilibrium
  delta   <- param["delta"]
  r       <- param["r"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  
  # compute roots of characteristic polynomials
  xi_1    <- no_commit_root2_fn(param)['xi_1']
  xi_2    <- no_commit_root2_fn(param)['xi_2']
  eta_1   <- commit_debt_root2_fn(param)['eta_1']
  eta_2   <- commit_debt_root2_fn(param)['eta_2']
  x_bar   <- constrained_reflecting_cutoff_fn(param,x_star)
  rho     <- x_star/x_bar
  
  # compute risk-free debt prices
  d_rf_delta <- (kappa+m)/(delta+m)
  d_rf_r     <- (kappa+m)/(r+m)
  
  # constants of integration in constrained region
  d_1     <- d_rf_r*eta_2*rho^eta_2/(eta_1*rho^eta_1 - eta_2*rho^eta_2)
  d_2     <- -d_rf_r*eta_1*rho^eta_1/(eta_1*rho^eta_1 - eta_2*rho^eta_2)
  
  # constants of integration in constrained region
  v_1     <- x_bar * (d_rf_delta*(1-xi_2*rho^(xi_2-1)) + 
                        xi_2*rho^(xi_2-1) / (delta-mu)/x_bar - 
                        d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
    (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
  v_2     <- x_bar * (d_rf_delta*(xi_1*rho^(xi_1-1) - 1) - 
                        xi_1*rho^(xi_1-1) / (delta-mu)/x_bar + 
                        d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
    (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
  
  # compute debt price when constrained
  debt  <- d_rf_r + d_1*(x/x_bar)^eta_1 + d_2*(x/x_bar)^eta_2
  
  # compute potential x_hat
  x_hat         <- xi_2*x_star/(1-xi_2)*
    ((1-xi_1)*v_1*rho^xi_1 + (1-xi_2)*v_2*rho^xi_2)/
    (xi_1*v_1*rho^xi_1 + xi_2*v_2*rho^xi_2)
  
  # compute debt price at x = x*
  debt_x_star  <- d_rf_r + d_1*rho^eta_1 + d_2*rho^eta_2
  value_x_star <- 1/(delta-mu)-d_rf_delta*x_star + v_1*rho^xi_1 + v_2*rho^xi_2
  value_x_hat  <- value_x_star + (x_star-x_hat) * debt_x_star
  
  # create flag indicating whether debt_x_star <> d_rf_delta
  region_3_flag     <- (debt_x_star < d_rf_delta)&(x_hat>0)&(x_hat<x_star)
  region_2_flag     <- 1-region_3_flag
  
  # compute debt price when x > x*
  ratio_1       <- x/x_bar
  d_high        <- d_rf_r + d_1*ratio_1^eta_1 + d_2*ratio_1^eta_2
  
  # compute debt price when x_low > x > x*
  d_mid         <- debt_x_star
  
  # compute debt price when x < x_low
  if (x_hat > 0){
    ratio_2       <- x/x_hat
    d_low         <- d_rf_delta - xi_2/x_hat*(ratio_2^(xi_2-1))*(value_x_hat+d_rf_delta*x_hat-1/(delta-mu))
  } else {
    d_low         <- 0
  }
  
  # compute debt value
  debt  <- (d_high*(x <= x_bar)*(x > x_star) + d_mid*(x > x_hat)*(x <= x_star) + 
              d_low*(x <= x_hat) + 0*(x > x_bar))*region_3_flag + 
    (d_high*(x <= x_bar)*(x > x_star) + d_mid*(x <= x_star) + 0*(x > x_bar))*region_2_flag
  
  return(debt)
}

constrained_reflecting_value_fn <- function(param, x_star, x){
  # compute value function for case with constrained for x > bar_x*
  # assume reflecting equilibrium
  # if d(x*) < d_rf_delta, 3 regions, including a smooth region for x in (0, hat_x)
  # if d(x*) > d_rf_delta, 2 regions only
  delta   <- param["delta"]
  r       <- param["r"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  
  # compute roots of characteristic polynomials
  xi_1    <- no_commit_root2_fn(param)['xi_1']
  xi_2    <- no_commit_root2_fn(param)['xi_2']
  eta_1   <- commit_debt_root2_fn(param)['eta_1']
  eta_2   <- commit_debt_root2_fn(param)['eta_2']
  x_bar   <- constrained_reflecting_cutoff_fn(param,x_star)
  rho     <- x_star/x_bar
  
  # compute risk-free debt prices
  d_rf_delta <- (kappa+m)/(delta+m)
  d_rf_r     <- (kappa+m)/(r+m)
  
  # constants of integration in constrained region
  v_1     <- x_bar * (d_rf_delta*(1-xi_2*rho^(xi_2-1)) + 
                        xi_2*rho^(xi_2-1) / (delta-mu)/x_bar - 
                        d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
    (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
  v_2     <- x_bar * (d_rf_delta*(xi_1*rho^(xi_1-1) - 1) - 
                        xi_1*rho^(xi_1-1) / (delta-mu)/x_bar + 
                        d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
    (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
  
  # compute value function when constrained
  value  <- 1/(delta-mu)-d_rf_delta*x + v_1*(x/x_bar)^xi_1 + v_2*(x/x_bar)^xi_2
  
  # constants of integration in constrained region for d(x)
  d_1     <- d_rf_r*eta_2*rho^eta_2/(eta_1*rho^eta_1 - eta_2*rho^eta_2)
  d_2     <- -d_rf_r*eta_1*rho^eta_1/(eta_1*rho^eta_1 - eta_2*rho^eta_2)
  
  # compute debt price and value function at x = x*
  debt_x_star   <- d_rf_r + d_1*rho^eta_1 + d_2*rho^eta_2
  value_x_star  <- 1/(delta-mu)-d_rf_delta*x_star + v_1*rho^xi_1 + v_2*rho^xi_2
  
  # compute potential x_hat
  x_hat         <- xi_2*x_star/(1-xi_2)*
    ((1-xi_1)*v_1*rho^xi_1 + (1-xi_2)*v_2*rho^xi_2)/
    (xi_1*v_1*rho^xi_1 + xi_2*v_2*rho^xi_2)
  
  # create flag indicating whether debt_x_star <> d_rf_delta
  region_3_flag     <- (debt_x_star < d_rf_delta)&(x_hat > 0)&(x_hat<x_star)
  region_2_flag     <- 1-region_3_flag
  
  # compute value at x_hat
  value_x_hat   <- value_x_star+(x_star-x_hat)*debt_x_star
  
  # compute value function when x > x*
  ratio_1       <- x/x_bar
  v_high        <- 1/(delta-mu) - d_rf_delta*x + v_1*ratio_1^xi_1+v_2*ratio_1^xi_2
  
  # compute value function when x_low > x > x*
  v_mid         <- value_x_star+(x_star-x)*debt_x_star
  
  # compute value function when x < x_hat
  if (x_hat > 0){
    ratio_2       <- x/x_hat
    v_low         <- 1/(delta-mu)*(1 - ratio_2^xi_2) - d_rf_delta*x*(1-ratio_2^(xi_2-1)) + value_x_hat*ratio_2^xi_2
  } else {
    v_low         <- 0
  }
  
  # return value function depending on case
  value  <- (v_high*(x <= x_bar)*(x > x_star) + v_mid*(x > x_hat)*(x <= x_star) + 
               v_low*(x <= x_hat) + 0*(x > x_bar))*region_3_flag + 
    (v_high*(x <= x_bar)*(x > x_star) + v_mid*(x <= x_star) + 0*(x > x_bar))*region_2_flag
  
  return(value)
}

constrained_hat_value_fn <- function(param, x_star, num_param){
  # compute value function for citizens in debt ceiling equilibrium
  # need to specify x_star, such that government cannot issue debt when x > x_star
  
  # load numerical parameters
  n_x       <- num_param['n_x']
  dt        <- num_param['dt']
  max_iter  <- num_param['max_iter']
  tol       <- num_param['tol']
  
  # load parameters
  delta     <- param["delta"]
  hat_delta <- param["hat_delta"]
  kappa     <- param["kappa"]
  m         <- param["m"]
  mu        <- param["mu"]
  sigma     <- param["sigma"]
  alpha     <- param["alpha"]
  theta     <- param["theta"]
  r         <- param["r"]
  nu        <- param["nu"]
  
  # set x_min so that we do not need to deal with singularity
  x_min    <- .01
  
  # check if equilibrium is reflecting or not
  x_cutoff_no_commit <- no_commit_cutoff_fn(param)
  x_star_min         <- issuance_constraint_ratio_limit_fn(param)*x_cutoff_no_commit
  
  if (x_star < x_star_min){
    # case with reflecting equilibrium
    x_bar <- constrained_reflecting_cutoff_fn(param, x_star)
    
    # compute roots with appropriate delta
    hat_param          <- param
    hat_param['delta'] <- param['hat_delta']
    xi_1               <- no_commit_root2_fn(hat_param)['xi_1']
    xi_2               <- no_commit_root2_fn(hat_param)['xi_2']
    eta_1              <- commit_debt_root2_fn(param)['eta_1']
    eta_2              <- commit_debt_root2_fn(param)['eta_2']
    xi_1b              <- no_commit_root2_fn(param)['xi_1']
    xi_2b              <- no_commit_root2_fn(param)['xi_2']
    rho                <- x_star/x_bar
    
    # compute risk-free debt prices
    d_rf_hat_delta <- (kappa+m)/(hat_delta+m)
    d_rf_delta     <- (kappa+m)/(delta+m)
    d_rf_r         <- (kappa+m)/(r+m)
    
    # constants of integration in constrained region, with hat_delta as discount rate
    v_1     <- x_bar * (d_rf_hat_delta*(1-xi_2*rho^(xi_2-1)) + 
                          xi_2*rho^(xi_2-1) / (hat_delta-mu)/x_bar - 
                          d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
      (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
    v_2     <- x_bar * (d_rf_hat_delta*(xi_1*rho^(xi_1-1) - 1) - 
                          xi_1*rho^(xi_1-1) / (hat_delta-mu)/x_bar + 
                          d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
      (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
    
    # constants of integration in constrained region for d(x)
    d_1     <- d_rf_r*eta_2*rho^eta_2/(eta_1*rho^eta_1 - eta_2*rho^eta_2)
    d_2     <- -d_rf_r*eta_1*rho^eta_1/(eta_1*rho^eta_1 - eta_2*rho^eta_2)
    
    # create grid and initial guess
    x     <- seq(from = x_min, to = x_bar, length.out = n_x)
    v     <- array(NA,n_x)
    
    # compute value at x_star, as well as debt price
    v_x_star  <- 1/(hat_delta-mu)-d_rf_hat_delta*x_star + v_1*rho^xi_1 + v_2*rho^xi_2
    d_x_star  <- d_rf_r + d_1*rho^eta_1 + d_2*rho^eta_2
    
    # compute value function when no issuance is allowed
    v[x > x_star]  <- 1/(hat_delta-mu) - d_rf_hat_delta*x[x > x_star] + 
      v_1*(x[x > x_star]/x_bar)^xi_1 + v_2*(x[x > x_star]/x_bar)^xi_2
    
    # determine whether there is (or not) a smooth region
    # first, compute constants of integration in constrained region using "correct" roots
    v_1b <- x_bar * (d_rf_delta*(1-xi_2b*rho^(xi_2b-1)) + 
                       xi_2b*rho^(xi_2b-1) / (delta-mu)/x_bar - 
                       d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
      (xi_1b*rho^(xi_1b-1) - xi_2b*rho^(xi_2b-1))
    v_2b <- x_bar * (d_rf_delta*(xi_1b*rho^(xi_1b-1) - 1) - 
                       xi_1b*rho^(xi_1b-1) / (delta-mu)/x_bar + 
                       d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
      (xi_1b*rho^(xi_1b-1) - xi_2b*rho^(xi_2b-1))
    
    # then, compute x_hat, "jump to" point
    x_hat <- xi_2b*x_star/(1-xi_2b)*
      ((1-xi_1b)*v_1b*rho^xi_1b + (1-xi_2b)*v_2b*rho^xi_2b)/
      (xi_1b*v_1b*rho^xi_1b + xi_2b*v_2b*rho^xi_2b)
    
    # check whether there is a smooth region near x = 0 or not -- depends on d_x_star
    if ((d_x_star < d_rf_delta)&(x_hat > 0)&(x_hat<x_star)){
      # case where there is a smooth region
      # compute value function in jump region
      v[(x <= x_star)&(x >= x_hat)] <- v_x_star + (x_star - x[(x <= x_star)&(x >= x_hat)])*d_x_star
      
      # compute value function at x_hat
      hat_v_x_hat <- v_x_star + (x_star-x_hat)*d_x_star
      
      # construct government value function in smooth region [0, x_hat], as well as financing policy
      value_x_star  <- 1/(delta-mu)-d_rf_delta*x_star + v_1b*rho^xi_1b + v_2b*rho^xi_2b
      value_x_hat   <- value_x_star+(x_star-x_hat)*d_x_star
      ratio_2       <- x[x < x_hat]/x_hat
      gov_v_p_low   <- -(xi_2b/x_hat)/(delta-mu)*ratio_2^(xi_2b-1) - 
        d_rf_delta*(1-ratio_2^(xi_2b-1)) + (xi_2b-1)*d_rf_delta*ratio_2^(xi_2b-1) + 
        xi_2b*(value_x_hat/x_hat)*ratio_2^(xi_2b-1)
      gov_v_pp_low  <- -(xi_2b*(xi_2b-1)/x_hat^2)/(delta-mu)*ratio_2^(xi_2b-2) + 
        ((xi_2b-1)/x_hat)*d_rf_delta*ratio_2^(xi_2b-2) + (((xi_2b-1)^2)/x_hat)*d_rf_delta*ratio_2^(xi_2b-2) + 
        (xi_2b-1)*xi_2b*(value_x_hat/x_hat^2)*ratio_2^(xi_2b-2)
      
      # construct financing policy
      g_low         <- (r-delta)*gov_v_p_low/gov_v_pp_low - nu*sigma*x[x < x_hat]
      d_low         <- - gov_v_p_low
      
      # use finite difference to compute value function in the smooth region
      x_low         <- x[x < x_hat]
      n_low         <- length(x_low)
      v_n_low       <- no_commit_value_fn(param, x_low)
      cons_low      <- 1 - (kappa+m)*x_low + g_low*d_low
      
      # compute dx
      dx     <- x[2]-x[1]
      
      # compute drift rate (invariant across iterations)
      drift  <- g_low - (mu + m)*x_low
      
      # construct matrix A
      A      <- Matrix(0, nrow = n_low, ncol = n_low, sparse = TRUE)
      
      # diagonal elements for A
      locs     <- cbind(seq.int(n_low),seq.int(n_low))
      A[locs]  <- -(abs(drift)/dx + (sigma*x_low/dx)^2)
      
      # lower diagonal elements -- probability to go to lower x
      locs     <- cbind((2:n_low),(1:(n_low-1)))
      A[locs]  <- -pmin(0,drift[-1])/dx + .5*(sigma*x_low[-1]/dx)^2
      
      # upper diagonal elements -- probability to go to higher x
      locs     <- cbind((1:(n_low-1)),(2:n_low))
      A[locs]  <- pmax(0,drift[-n_low])/dx + .5*(sigma*x_low[-n_low]/dx)^2
      
      # impose reflection at x = 0 and x = x_bar
      A[1,1]         <- -sum(A[1,-1])
      A[n_low,n_low] <- -sum(A[n_low,-n_low])
      
      # construct B matrix -- it assumes value matching at x = x_hat
      B                  <- (1/dt + hat_delta - mu)*Diagonal(n_low) - A
      B[n_low,n_low]     <- 1
      B[n_low,(n_low-1)] <- 0
      
      # initialization
      err     <- 1      # initialize error for value function
      n_iter  <- 1      # initialize iteration #
      while ((err>tol)&(n_iter<max_iter)){
        # if error greater than tolerance
        # assign most recent value to v
        v_low <- v_n_low
        
        # enter right hand side, including correction for boundary conditions
        RHS          <- cons_low + v_low/dt
        RHS[n_low]   <- hat_v_x_hat
        
        # solve linear system
        v_n_low  <- c(array(solve(B,RHS)))
        
        # compute error
        err      <- max(abs(v_n_low-v_low)/dt)
        
        # update iteration #
        n_iter  <- n_iter + 1
      }
      
      # compute value function in smooth region x < x_low
      v[x < x_hat] <- v_n_low
      
      # interpolate at zero
      hat_v_0  <- v[1] - x_min*(v[2]-v[1])/dx
      x        <- c(0,x)
      v        <- c(hat_v_0,v)
    } else {
      # case where there is no smooth region and there is only a jump from 0 to X*
      v[x <= x_star] <- v_x_star + (x_star - x[x <= x_star])*d_x_star
      # compute v_0
      hat_v_0  <- v_x_star + x_star*d_x_star
      x        <- c(0,x)
      v        <- c(hat_v_0,v)
    }
  } else {
    # case with smooth equilibrium
    x_bar <- no_commit_cutoff_fn(param)
    
    # create grid
    x     <- seq(from = x_min, to = x_bar, length.out = n_x)
    
    # compute debt price and issuance policy, as well as initial guess for hat v
    g               <- no_commit_g_fn(param, x)
    g[x>x_star]     <- 0
    cons            <- no_commit_c_fn(param, x)
    cons[x>x_star]  <- 1 - (kappa+m)*x[x>x_star]
    
    # guess value function
    orig_v <- no_commit_value_fn(param, x)
    v_n    <- orig_v
    
    # compute dx
    dx     <- x[2]-x[1]
    
    # compute drift rate (invariant across iterations)
    drift  <- g - (mu + m)*x
    
    # construct matrix A
    A      <- Matrix(0, nrow = n_x, ncol = n_x, sparse = TRUE)
    
    # diagonal elements for A
    locs   <- cbind(seq.int(n_x),seq.int(n_x))
    A[locs]  <- -(abs(drift)/dx + (sigma*x/dx)^2)
    
    # lower diagonal elements -- probability to go to lower x
    locs     <- cbind((2:n_x),(1:(n_x-1)))
    A[locs]  <- -pmin(0,drift[-1])/dx + .5*(sigma*x[-1]/dx)^2
    
    # upper diagonal elements -- probability to go to higher x
    locs     <- cbind((1:(n_x-1)),(2:n_x))
    A[locs]  <- pmax(0,drift[-n_x])/dx + .5*(sigma*x[-n_x]/dx)^2
    
    # impose reflection at x = 0 and x = x_bar
    A[1,1]     <- -sum(A[1,-1])
    A[n_x,n_x] <- -sum(A[n_x,-n_x])
    
    # construct B matrix
    B              <- (1/dt + hat_delta - mu)*Diagonal(n_x) - A
    
    # boundary condition at x = x_bar
    idx_reentry        <- which.min(abs(x-theta*x_bar))
    B[n_x,]            <- 0
    B[n_x,n_x]         <- 1
    B[n_x,idx_reentry] <- -alpha
    
    # constant part of the right-handside + take care of boundary conditions
    RHS_cst            <- cons
    
    # initialization
    err     <- 1      # initialize error for value function
    err_vec <- c(1)
    iter    <- 1      # initialize iteration #
    # compute citizens' welfare value
    while ((err>tol)&(iter<max_iter)){
      # if error greater than tolerance
      # assign most recent value to v
      v   <- v_n
      
      # enter right hand side, including correction for boundary conditions
      RHS      <- RHS_cst + v/dt
      RHS[n_x] <- 0
      
      # solve linear system
      v_n  <- c(array(solve(B,RHS)))
      
      # compute error
      err      <- max(abs(v_n-v)/dt)
      err_vec  <- c(err_vec,err)
      
      # update iteration #
      iter     <- iter + 1
    }
    
    # interpolate at zero
    hat_v_0  <- v_n[1] - x_min*(v_n[2]-v_n[1])/dx
    x        <- c(0,x)
    v        <- c(hat_v_0,v_n)
  }
  
  # smooth function
  hat_v_fn <- approxfun(x, v, method = "linear")
  
  # output
  output        <- list(hat_v_fn,hat_v_0)
  names(output) <- c("hat_v", "hat_v0")
  return(output)
}

constrained_reflecting_cutoff_fn <- function(param, x_star){
  # compute optimal default cutoff for constrained reflecting equilibrium
  delta   <- param["delta"]
  r       <- param["r"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  
  # compute roots of characteristic polynomials
  xi_1    <- no_commit_root2_fn(param)['xi_1']
  xi_2    <- no_commit_root2_fn(param)['xi_2']
  eta_1   <- commit_debt_root2_fn(param)['eta_1']
  eta_2   <- commit_debt_root2_fn(param)['eta_2']
  
  # compute risk-free debt prices
  d_rf_delta <- (kappa+m)/(delta+m)
  d_rf_r     <- (kappa+m)/(r+m)
  
  x_bar_fn <- function(x_bar_test){
    # implicit function that will be used to compute x_bar
    rho   <- x_star/x_bar_test
    # constants of integration in constrained region
    v_1   <- x_bar_test * (d_rf_delta*(1-xi_2*rho^(xi_2-1)) +
                             xi_2*rho^(xi_2-1) / (delta-mu)/x_bar_test -
                             d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
      (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
    v_2   <- x_bar_test * (d_rf_delta*(xi_1*rho^(xi_1-1) - 1) - 
                             xi_1*rho^(xi_1-1) / (delta-mu)/x_bar_test +
                             d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
      (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
    # compute output
    output<- -d_rf_delta*x_bar_test + v_1 * xi_1 + v_2 * xi_2
    return(output)
  }
  
  # compute guess for x_bar
  x_bar_guess <- no_commit_cutoff_fn(param)
  
  # solve implicit equation
  x_bar       <- uniroot(x_bar_fn, lower = x_bar_guess, upper = 10*x_bar_guess, maxiter = 1000)$root
  
  return(x_bar)
}

mult_eq_fn <- function(param, x_star, x_bar_test){
  # compute optimal default cutoff for constrained reflecting equilibrium
  delta   <- param["delta"]
  r       <- param["r"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  
  # compute roots of characteristic polynomials
  xi_1    <- no_commit_root2_fn(param)['xi_1']
  xi_2    <- no_commit_root2_fn(param)['xi_2']
  eta_1   <- commit_debt_root2_fn(param)['eta_1']
  eta_2   <- commit_debt_root2_fn(param)['eta_2']
  
  # compute risk-free debt prices
  d_rf_delta <- (kappa+m)/(delta+m)
  d_rf_r     <- (kappa+m)/(r+m)
  
  # implicit function that will be used to compute x_bar
  rho   <- x_star/x_bar_test
  
  # constants of integration in constrained region
  v_1   <- x_bar_test * (d_rf_delta*(1-xi_2*rho^(xi_2-1)) +
                           xi_2*rho^(xi_2-1) / (delta-mu)/x_bar_test -
                           d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
    (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
  v_2   <- x_bar_test * (d_rf_delta*(xi_1*rho^(xi_1-1) - 1) - 
                           xi_1*rho^(xi_1-1) / (delta-mu)/x_bar_test +
                           d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
    (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
  
  # compute output
  output<- -d_rf_delta*x_bar_test + v_1 * xi_1 + v_2 * xi_2
  
  return(output)
}

mult_eq_fn2 <- function(param, x_star, x_bar_test){
  # compute optimal default cutoff for constrained reflecting equilibrium
  delta   <- param["delta"]
  r       <- param["r"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  
  # compute roots of characteristic polynomials
  xi_1    <- no_commit_root2_fn(param)['xi_1']
  xi_2    <- no_commit_root2_fn(param)['xi_2']
  eta_1   <- commit_debt_root2_fn(param)['eta_1']
  eta_2   <- commit_debt_root2_fn(param)['eta_2']
  
  # compute risk-free debt prices
  d_rf_delta <- (kappa+m)/(delta+m)
  d_rf_r     <- (kappa+m)/(r+m)
  
  # implicit function that will be used to compute x_bar
  rho   <- x_star/x_bar_test
  
  # constants of integration in constrained region
  denom <- rho^(xi_2-1) - rho^(xi_1-1)
  num1  <- (xi_2-xi_1)*(d_rf_delta-d_rf_r)
  num2  <- (xi_1*(1-xi_2)*rho^(xi_1-1)-
              xi_2*(1-xi_1)*rho^(xi_2-1))*d_rf_delta
  num3  <- -(xi_2-xi_1)*d_rf_r*(eta_2-eta_1)/
    (eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))
  num4  <- -xi_1*xi_2/(x_bar_test*(delta-mu))
  
  # compute output
  output<- (num1+num2+num3)/denom+num4
  
  return(output)
}

constrained_reflecting_jump_cutoff_fn <- function(param, x_star){
  # compute jump cutoff for constrained reflecting equilibrium
  delta   <- param["delta"]
  r       <- param["r"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  
  # compute roots of characteristic polynomials
  xi_1    <- no_commit_root2_fn(param)['xi_1']
  xi_2    <- no_commit_root2_fn(param)['xi_2']
  eta_1   <- commit_debt_root2_fn(param)['eta_1']
  eta_2   <- commit_debt_root2_fn(param)['eta_2']
  
  # compute risk-free debt prices
  d_rf_delta <- (kappa+m)/(delta+m)
  d_rf_r     <- (kappa+m)/(r+m)
  
  # compute default cutoff
  x_bar      <- constrained_reflecting_cutoff_fn(param, x_star)
  
  ushort_x_fn <- function(ushort_x_test){
    # implicit function that will be used to compute ushort_x
    rho   <- x_star/x_bar
    
    # constants of integration for value function in constrained region
    v_1   <- x_bar * (d_rf_delta*(1-xi_2*rho^(xi_2-1)) +
                        xi_2*rho^(xi_2-1) / (delta-mu)/x_bar -
                        d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
      (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
    v_2   <- x_bar * (d_rf_delta*(xi_1*rho^(xi_1-1) - 1) - 
                        xi_1*rho^(xi_1-1) / (delta-mu)/x_bar +
                        d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
      (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
    
    # constants of integration for debt price in constrained region
    d_1   <- d_rf_r*eta_2*(rho^eta_2)/(eta_1*rho^eta_1 - eta_2*rho^eta_2)
    d_2   <- -d_rf_r*eta_1*(rho^eta_1)/(eta_1*rho^eta_1 - eta_2*rho^eta_2)
    
    # compute d(x_star), v(x_star) and v(ushort_x)
    d_star    <- d_rf_r + d_1*rho^eta_1 + d_2*rho^eta_2
    v_star    <- 1/(delta-mu) - d_rf_delta*x_star + v_1*rho^xi_1 + v_2*rho^xi_2
    v_ushort  <- v_star+(x_star - ushort_x_test)*d_star
    
    # compute d(x_star2)
    d_star2   <- d_rf_delta - xi_2*(v_ushort/ushort_x_test + d_rf_delta - 1/(delta-mu)/ushort_x_test)
    
    # compute output
    output<- d_star-d_star2
    return(output)
  }
  
  ushort_x_guess <- .0001
  ushort_x       <- uniroot(ushort_x_fn, lower = ushort_x_guess, upper = x_star, maxiter = 1000)$root
  
  return(ushort_x)
}

constrained_reflecting_jump_cutoff_fn2 <- function(param, x_star){
  # compute jump cutoff for constrained reflecting equilibrium
  # algebra intensive strategy
  delta   <- param["delta"]
  r       <- param["r"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  
  # compute roots of characteristic polynomials
  xi_1    <- no_commit_root2_fn(param)['xi_1']
  xi_2    <- no_commit_root2_fn(param)['xi_2']
  eta_1   <- commit_debt_root2_fn(param)['eta_1']
  eta_2   <- commit_debt_root2_fn(param)['eta_2']
  
  # compute risk-free debt prices
  d_rf_delta <- (kappa+m)/(delta+m)
  d_rf_r     <- (kappa+m)/(r+m)
  
  # compute default cutoff
  x_bar      <- constrained_reflecting_cutoff_fn(param, x_star)
  
  # implicit function that will be used to compute ushort_x
  rho   <- x_star/x_bar
  
  # constants of integration for value function in constrained region
  v_1   <- x_bar * (d_rf_delta*(1-xi_2*rho^(xi_2-1)) +
                      xi_2*rho^(xi_2-1) / (delta-mu)/x_bar -
                      d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
    (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
  v_2   <- x_bar * (d_rf_delta*(xi_1*rho^(xi_1-1) - 1) - 
                      xi_1*rho^(xi_1-1) / (delta-mu)/x_bar +
                      d_rf_r*(1+(eta_2-eta_1)/(eta_1*rho^(-eta_2) - eta_2*rho^(-eta_1))))/
    (xi_1*rho^(xi_1-1) - xi_2*rho^(xi_2-1))
  
  # constants of integration for debt price in constrained region
  d_1   <- d_rf_r*eta_2*(rho^eta_2)/(eta_1*rho^eta_1 - eta_2*rho^eta_2)
  d_2   <- -d_rf_r*eta_1*(rho^eta_1)/(eta_1*rho^eta_1 - eta_2*rho^eta_2)
  
  # compute ushort x val
  ushort_x   <- xi_2*x_star/(1-xi_2)*
    ((1-xi_1)*v_1*rho^xi_1 + (1-xi_2)*v_2*rho^xi_2)/
    (xi_1*v_1*rho^xi_1 + xi_2*v_2*rho^xi_2)
  
  return(ushort_x)
}

#===================================
# Issuance Cap Equilibrium Functions
#===================================

max_g_eq_fn <- function(param, num_param){
  # compute optimal boundaries
  # solve for v() continuous at x = x* and g() continuous at x = x*
  # for given bar x:         ## if min g < 0, increase x* 
  ## if g jump < 0, increase x* 
  ## if g jump > 0, decrease x* 
  # if v jump > 0, increase bar x
  # if v jump < 0, increase bar x
  
  # specify maximum number of possible tests of equilibrium
  max_iter_v     <- 20
  max_iter_g     <- 20
  
  # load parameters
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  bar_g   <- param["bar_g"]
  
  # compute x_bar when there is no commitment
  x_bar_no_commit <- no_commit_cutoff_fn(param)
  
  # error tolerance
  v_err_tol  <- .01
  g_err_tol  <- .01
  
  # initialize errors + tolerance
  v_error       <- 10
  v_error_vec   <- v_error
  v_iter        <- 1
  while ((abs(v_error) > v_err_tol)&(v_iter < max_iter_v)){
    # update x_bar_guess
    if (v_iter == 1) {
      x_bar_guess      <- x_bar_no_commit
      x_bar_guess_vec  <- c(0,x_bar_no_commit)
    } else if ((prod(v_error_vec>0) == 1)&(v_iter>1)) {
      x_bar_guess      <- 2*x_bar_guess_vec[length(x_bar_guess_vec)]
      x_bar_guess_vec  <- c(x_bar_guess_vec,x_bar_guess)
    } else if ((prod(v_error_vec>0) == 0)&(v_iter>1)) {
      step_size        <- .5*abs(x_bar_guess_vec[(length(x_bar_guess_vec)-1)] - x_bar_guess_vec[length(x_bar_guess_vec)])
      x_bar_guess      <- x_bar_guess_vec[length(x_bar_guess_vec)] + step_size*(1*(v_error > 0) - 1*(v_error < 0))
      x_bar_guess_vec  <- c(x_bar_guess_vec,x_bar_guess)
    }
    
    # begin loop
    g_error     <- 10
    g_iter      <- 1
    g_error_vec <- 1
    while ((abs(g_error) > g_err_tol)&(g_iter < max_iter_g)){
      # update mult_guess
      if (g_iter == 1) {
        mult_guess    <- .5
        mult_vec      <- .5
        x_star_guess  <- mult_guess*x_bar_guess
      } else {
        mult_guess    <- mult_vec[length(mult_vec)] + (0.5^g_iter)*(1*((g_error < 0)|(min_g < -.01)) - 
                                                                      1*((g_error > 0)&(min_g > -.01)))
        x_star_guess  <- mult_guess*x_bar_guess
        mult_vec      <- c(mult_vec, mult_guess)
      }
      
      # recompute equilibrium
      dummy_res   <- max_g_partial_eq_fn(param, x_star = x_star_guess, x_bar = x_bar_guess, num_param)
      g_error     <- dummy_res[['g_jump']]
      g_error_vec <- c(g_error_vec,g_error)
      min_g       <- dummy_res[['min_g']]
      g_iter      <- g_iter + 1
    }
    
    # recompute equilibrium
    v_error      <- dummy_res[['v_jump']]
    v_error_vec  <- c(v_error_vec, v_error)
    v_iter       <- v_iter + 1
  }
  
  eq_output <- list(dummy_res[['v']], dummy_res[['d']], 
                    dummy_res[['g']], dummy_res[['c']], dummy_res[['incentive']],
                    x_star_guess,x_bar_guess, v_error, g_error)
  
  names(eq_output) <- c("v", "d", "g", "c", "incentive", "x_star", "x_bar", "error_v", "error_g")
  return(eq_output)
}

max_g_partial_eq_fn <- function(param, x_star, x_bar, num_param){
  # compute debt price and value function for case without commitment
  # if speed_flag == 1, plot issuance policy, consumption policy, value function and debt price
  # maximum issuance rate constraint
  # compute for given assumed boundaries
  # we impose (a) debt price is continuous, (b) value function is zero at default boundary, (c) smooth pasting at default boundary
  # what is missing from equilibrium: (a) continuity of v, and (b) continuity of issuance policy
  
  # folder dump
  dump_dir    <- "~/Dropbox/Research/Sovereign Defaults and Commitment/Code/Dump"
  speed_flag  <- 1
  
  # numerical solution parameters
  n_x        <- num_param['n_x']
  dt         <- num_param['dt']
  max_iter   <- num_param['max_iter']
  tol        <- num_param['tol']
  
  # load parameters
  delta   <- param["delta"]
  kappa   <- param["kappa"]
  m       <- param["m"]
  mu      <- param["mu"]
  sigma   <- param["sigma"]
  alpha   <- param["alpha"]
  theta   <- param["theta"]
  r       <- param["r"]
  bar_g   <- param["bar_g"]
  
  # compute risk-free debt price
  d_rf     <- (kappa+m)/(delta+m)
  
  # compute root of characteristic polynomial
  xi      <- no_commit_root2_fn(param)
  xi_1    <- xi["xi_1"]
  xi_2    <- xi["xi_2"]
  
  # compute the coefficients v_1 and v_2
  # those coefficients are computed such that there is value matching and smooth pasting at x = bar x
  v_1      <- (xi_2-1)/(xi_2-xi_1)*d_rf*x_bar - xi_2/(xi_2-xi_1)/(delta-mu)
  v_2      <- (1-xi_1)/(xi_2-xi_1)*d_rf*x_bar + xi_1/(xi_2-xi_1)/(delta-mu)
  
  # compute the coefficients d_1 and d_2
  d_1      <- -v_1/x_bar*xi_1
  d_2      <- -v_2/x_bar*xi_2
  
  # compute value function at boundary x^*
  v_x_star       <- 1/(delta - mu) - x_star*d_rf + v_1*(x_star/x_bar)^xi_1 + v_2*(x_star/x_bar)^xi_2
  
  # compute debt price at boundary x^*
  d_x_star       <- d_rf + d_1*(x_star/x_bar)^(xi_1-1) + d_2*(x_star/x_bar)^(xi_2-1)
  
  # compute debt price between x = 0 and x = x_star
  x      <- seq(from = 0, to = x_star, length.out = n_x)
  dx     <- x[2]-x[1]
  d_next <- seq(from = (kappa+m)/(r+m), to = d_x_star, length.out = n_x)
  
  # compute drift rate (invariant across iterations)
  drift  <- bar_g - (mu + m - sigma^2)*x
  disc   <- rep(r+m,n_x)
  
  # construct diagonals of matrix to invert for debt computations
  A_diag     <- 1 + dt*(disc+abs(drift)/dx+(sigma*x/dx)^2)
  A_up_diag  <- -dt*(pmax(drift,0)/dx+0.5*(sigma*x/dx)^2)
  A_lo_diag  <- -dt*(-pmin(drift,0)/dx+0.5*(sigma*x/dx)^2)
  
  # construct matrix to invert for debt prices
  A_mat                                <- Matrix(0, nrow = n_x, ncol = n_x, sparse = TRUE)
  A_mat[cbind((1:(n_x-1)),(2:(n_x)))]  <- A_up_diag[-n_x]
  A_mat[cbind((2:(n_x)),(1:(n_x-1)))]  <- A_lo_diag[-1]
  A_mat[cbind((1:n_x),(1:n_x))]        <- A_diag
  
  # boundary conditions for debt price -- natural boundary at 0 and value matching at x = x_star
  A_mat[n_x,(n_x-1)]  <- 0
  A_mat[n_x,n_x]      <- 1
  
  # compute the constant part of the RHS
  RHS_cst             <- c(array(dt*(kappa+m), n_x))
  
  # initialization
  err_d     <- 1      # initialize error for debt
  err_v     <- 1      # initialize error for value function
  n_iter    <- 1      # initialize iteration #
  # compute debt price
  while ((err_d>tol)&(n_iter<max_iter)){
    # if error greater than tolerance
    # assign most recent value to d
    d   <- d_next
    
    # enter right-handside, including correction for boundary at x = x_star
    RHS      <- RHS_cst+d
    RHS[n_x] <- d_x_star
    
    # solve linear system
    d_next  <- solve(A_mat,RHS)
    d_next  <- c(array(d_next))
    
    # compute error
    err_d <- max(abs(d_next-d)/dt)
    
    # update iteration #
    n_iter  <- n_iter + 1
  }
  
  # update debt price
  d     <- d_next
  
  # compute value function
  v_next    <- seq(from = 1/(delta - mu), to = v_x_star, length.out = n_x)
  
  # reinitialize iteration count
  n_iter    <- 1
  
  # compute drift rate (invariant across iterations) and utility flow rate
  drift <- bar_g - (mu + m)*x
  disc  <- rep(delta-mu,n_x)
  flow  <- 1 + bar_g * d - (kappa+m)*x
  
  # construct diagonals of matrix to invert for value function computations
  B_up_diag  <- -dt*(pmax(drift,0)/dx+0.5*(sigma*x/dx)^2)
  B_lo_diag  <- -dt*(-pmin(drift,0)/dx+0.5*(sigma*x/dx)^2)
  B_diag     <- 1 + dt*(disc+abs(drift)/dx+(sigma*x/dx)^2)
  
  # construct matrix to invert for value function prices
  B_mat                                <- Matrix(0, nrow = n_x, ncol = n_x, sparse = TRUE)
  B_mat[cbind((1:(n_x-1)),(2:(n_x)))]  <- B_up_diag[-n_x]
  B_mat[cbind((2:(n_x)),(1:(n_x-1)))]  <- B_lo_diag[-1]
  B_mat[cbind((1:(n_x)),(1:(n_x)))]    <- B_diag
  
  # boundary conditions for debt price -- natural boundary at 0 and value matching at x = x_star
  B_mat[n_x,(n_x-1)]  <- -1/dx
  B_mat[n_x,n_x]      <- 1/dx
  
  # boundary conditions -- impose smooth pasting at x = x_star
  RHS[n_x] <- -d_x_star
  
  # compute value function
  while ((err_v>tol)&(n_iter<max_iter)){
    # if error greater than tolerance
    # assign most recent value to v
    v            <- v_next
    
    # enter right-handside, including correction for boundary at x = x_star
    RHS[-n_x] <- dt*flow[-n_x]+v[-n_x]
    
    # solve linear system
    v_next  <- solve(B_mat,RHS)
    v_next  <- c(array(v_next))
    
    # compute error
    err_v <- max(abs(v_next-v)/dt)
    
    # update iteration #
    n_iter  <- n_iter + 1
  }
  
  # compute debt price and value function across the state space
  x_add <- seq(from = x_star+dx, to = x_bar, length.out = n_x)
  v_add <- 1/(delta - mu)-d_rf*x_add + v_1*(x_add/x_bar)^xi_1 + v_2*(x_add/x_bar)^xi_2
  d_add <- d_rf + d_1*(x_add/x_bar)^(xi_1-1) + d_2*(x_add/x_bar)^(xi_2-1)
  
  # compute incentive d(x)+v'(x) in constrained region
  v_prime   <- diff(v)/dx
  v_prime   <- c(v_prime[1], v_prime)
  inc       <- d + v_prime
  inc_add   <- rep(0,n_x)
  
  # compute issuance and consumption policy in constrained region
  g    <- rep(bar_g,n_x)
  c    <- 1 + bar_g*d - (kappa+m)*x
  
  # compute issuance and consumption policy in unconstrained region
  d_prime_add <- (xi_1-1)*(d_1/x_bar)*(x_add/x_bar)^(xi_1-2) + (xi_2-1)*(d_2/x_bar)*(x_add/x_bar)^(xi_2-2)
  g_add       <- (r - delta)*d_add/d_prime_add
  c_add       <-  1 - (kappa+m)*x_add + g_add*d_add
  
  # compute discontinuity in value function and in issuance policy
  v_jump   <- v[n_x] - v_add[1]
  g_jump   <- bar_g - g_add[1]
  min_g    <- min(g_add)
  
  if (speed_flag == 0){
    # if we need to dump the plots
    setwd(dump_dir)
    x_bar_short  <- signif(x_bar, digits = 3)
    x_star_short <- signif(x_star, digits = 3)
    
    # plot value function
    file_name <- paste("value_x_star_=_",x_star_short,"_x_bar_=_",x_bar_short,".png",sep="")
    png(paste(file_name))
    plot(x = c(x,x_add), y = c(v_next,v_add), bty = 'n', type = 'l', xlab = 'x', ylab = 'v')
    abline(h = 1/(delta-mu), col = 'green', lty = 2)
    abline(v = x_star, col = 'red', lty = 2)
    dev.off()
    
    # plot debt price
    file_name <- paste("debt_x_star_=_",x_star_short,"_x_bar_=_",x_bar_short,".png",sep="")
    png(paste(file_name))
    plot(x = c(x,x_add), y = c(d_next,d_add), bty = 'n', type = 'l', xlab = 'x', ylab = 'd')
    abline(v = x_star, col = 'red', lty = 2)
    dev.off()
    
    # plot issuance policy
    file_name <- paste("g_x_star_=_",x_star_short,"_x_bar_=_",x_bar_short,".png",sep="")
    png(paste(file_name))
    plot(x = c(x,x_add), y = c(g,g_add), bty = 'n', type = 'l', xlab = 'x', ylab = 'g')
    abline(v = x_star, col = 'red', lty = 2)
    dev.off()
    
    # plot consumption policy
    file_name <- paste("cons_x_star_=_",x_star_short,"_x_bar_=_",x_bar_short,".png",sep="")
    png(paste(file_name))
    plot(x = c(x,x_add), y = c(c,c_add), bty = 'n', type = 'l', xlab = 'x', ylab = 'consumption')
    abline(v = x_star, col = 'red', lty = 2)
    dev.off()
    
    # plot incentives
    file_name <- paste("incentives_x_star_=_",x_star_short,"_x_bar_=_",x_bar_short,".png",sep="")
    png(paste(file_name))
    plot(x = c(x,x_add), y = c(inc,inc_add), bty = 'n', type = 'l', xlab = 'x', ylab = 'incentives')
    abline(v = x_star, col = 'red', lty = 2)
    dev.off()
  }
  
  # create value function and debt price function
  x    <- c(x,x_add)
  v    <- c(v_next,v_add)
  d    <- c(d_next,d_add)
  g    <- c(g,g_add)
  c    <- c(c,c_add)
  inc  <- c(inc,inc_add)
  
  # smooth function
  d_fn     <- approxfun(x, d, method = "linear")
  v_fn     <- approxfun(x, v, method = "linear")
  g_fn     <- approxfun(x, g, method = "linear")
  c_fn     <- approxfun(x, c, method = "linear")
  inc_fn   <- approxfun(x, inc, method = "linear")
  
  # output
  output        <- list(v_fn, d_fn, g_fn, c_fn, inc_fn, v_jump, g_jump, min_g)
  names(output) <- c("v", "d", "g", "c", "incentive", "v_jump", "g_jump", "min_g")
  return(output)
}

commit_hat_value_fn <- function(param, user_spec, x, g, cons, num_param){
  # compute value function for citizens
  # compute for given assumed boundaries
  # if user_spec == 1, use x, g, cons for consumption and issuance policy
  # if user_spec == 0, recompute equilibrium
  
  # load numerical parameters
  n_x       <- num_param['n_x']
  dt        <- num_param['dt']
  max_iter  <- num_param['max_iter']
  tol       <- num_param['tol']
  
  # load parameters
  delta     <- param["delta"]
  hat_delta <- param["hat_delta"]
  kappa     <- param["kappa"]
  m         <- param["m"]
  mu        <- param["mu"]
  sigma     <- param["sigma"]
  alpha     <- param["alpha"]
  theta     <- param["theta"]
  r         <- param["r"]
  
  if (user_spec == 0) {
    # compute constrained equilibrium solution
    const_eq <- max_g_eq_fn(param)
    # compute default boundary without commitment
    x_bar    <- const_eq[['x_bar']]
    # set x_min so that we do not need to deal with singularity
    x_min    <- .01
    # compute debt price between x = 0 and x = x_star
    x        <- seq(from = x_min, to = x_bar, length.out = n_x)
    # compute debt price and issuance policy, as well as initial guess for hat v
    g        <- const_eq[['g']](x)
    cons     <- const_eq[['c']](x)
  } else {
    # just specify the default boundary
    x_bar    <- max(x)
    x_min    <- min(x)
    n_x      <- length(x)
  }
  
  # compute dx
  dx        <- x[2] - x[1]
  
  # compute initial guess
  orig_v    <- no_commit_value_fn(param, x)
  v_n       <- orig_v
  
  # compute drift rate (invariant across iterations)
  drift <- g - (mu + m)*x
  
  # construct matrix A
  A        <- Matrix(0, nrow = n_x, ncol = n_x, sparse = TRUE)
  
  # diagonal elements for A
  locs     <- cbind(seq.int(n_x),seq.int(n_x))
  A[locs]  <- -(abs(drift)/dx + (sigma*x/dx)^2)
  
  # lower diagonal elements -- probability to go to lower x
  locs     <- cbind((2:n_x),(1:(n_x-1)))
  A[locs]  <- -pmin(0,drift[-1])/dx + .5*(sigma*x[-1]/dx)^2
  
  # upper diagonal elements -- probability to go to higher x
  locs     <- cbind((1:(n_x-1)),(2:n_x))
  A[locs]  <- pmax(0,drift[-n_x])/dx + .5*(sigma*x[-n_x]/dx)^2
  
  # impose reflection at x = 0 and x = x_bar
  A[1,1]     <- -sum(A[1,-1])
  A[n_x,n_x] <- -sum(A[n_x,-n_x])
  
  # construct B matrix
  B                  <- (1/dt + hat_delta - mu)*Diagonal(n_x) - A
  idx_reentry        <- which.min(abs(x-theta*x_bar))
  B[n_x,n_x]         <- 1
  B[n_x,(n_x-1)]     <- 0
  B[n_x,idx_reentry] <- -alpha
  
  # constant part of the right-handside + take care of boundary conditions
  RHS_cst   <- cons
  
  # initialization
  err     <- 1      # initialize error for value function
  err_vec <- c(1)
  n_iter  <- 1      # initialize iteration #
  # compute debt price
  while ((err>tol)&(n_iter<max_iter)){
    # if error greater than tolerance
    # assign most recent value to d_vec
    v   <- v_n
    
    # enter right-handside, including correction for boundary conditions
    RHS      <- RHS_cst+v/dt
    RHS[n_x] <- 0
    
    # solve linear system
    v_n  <- c(array(solve(B,RHS)))
    
    # compute error
    err      <- max(abs(v_n-v)/dt)
    err_vec  <- c(err_vec,err)
    
    # update iteration #
    n_iter  <- n_iter + 1
  }
  
  # interpolate at zero
  hat_v_0  <- v_n[1] - x_min*(v_n[2]-v_n[1])/dx
  x        <- c(0,x)
  v        <- c(hat_v_0,v_n)
  
  # smooth function
  hat_v_fn <- approxfun(x, v, method = "linear")
  
  # output
  output        <- list(hat_v_fn,hat_v_0)
  names(output) <- c("hat_v", "hat_v0")
  return(output)
}

#================================
# Aguiar & Amador Model Functions
#================================

# parametrize density of outside option
g_out  <- function(x, param){
  V_high <- param['V_high']
  beta   <- param['beta']
  return(exp(-abs(x-V_high)/beta)/(2*beta))
}

# parametrize CDF of outside option
G_out  <- function(x, param){
  V_high <- param['V_high']
  beta   <- param['beta']
  out_1 <- exp((x - V_high)/beta)/2
  out_2 <- 1 - exp(-(x - V_high)/beta)/2
  return((x < V_high)*out_1 + (x >= V_high)*out_2)
}

# Smooth MPE in Amador & Aguiar Model
no_commit_AA_fn <- function(param, num_param){
  # model such that option to default is always V_low
  # with Poisson arrival time, default value is V_tilde, random variable with mean V_high
  
  # load parameters
  V_high <- param['V_high']
  V_low  <- param['V_low']
  y      <- param['y']
  beta   <- param['beta']
  delta  <- param['delta']
  lambda <- param['lambda']
  r      <- param['r']
  m      <- param['m']
  kappa  <- param['kappa']
  
  # load numerical parameters
  n_x    <- num_param['n_x']
  n_F    <- num_param['n_F']
  F_max  <- num_param['F_max']
  F_min  <- num_param['F_min']
  
  # numerical parameters
  dt       <- num_param['dt']
  max_iter <- num_param['max_iter']
  tol      <- num_param['tol']
  conv_min <- num_param['conv_min']
  
  # create x grid for random value 
  x_min  <- V_high - 10*sqrt(2)*beta
  x_max  <- V_high + 10*sqrt(2)*beta
  x_grid <- seq(from = x_min, to = x_max, length.out = n_x)
  dx     <- x_grid[2]-x_grid[1]
  g_vec  <- g_out(x_grid, param)
  
  # create F_grid
  F_vec <- seq(from = F_min, to = F_max, length.out = n_F)
  dF    <- F_vec[2]-F_vec[1]
  
  # compute asymptotic value of V as F = F_max
  D_min  <- (kappa + m)/(delta + lambda + m)
  D_max  <- (kappa + m)/(delta + m)
  V_asym <- (y + lambda * V_high) / (delta + lambda) - D_min*F_vec
  
  # construct transition intensity matrix
  A       <- Matrix(0, nrow = n_F, ncol = n_F, sparse = TRUE)
  # diagonal elements for A
  locs    <- cbind(seq.int(n_F),seq.int(n_F))
  A[locs] <- -m*F_vec/dF
  # lower diagonal elements -- probability to go to lower F
  locs     <- cbind((2:n_F),(1:(n_F-1)))
  A[locs]  <- -pmin(0,-m*F_vec[-1]/dF)
  # impose reflection at F = 0
  A[1,2]   <- -A[1,1]
  # check that matrix is correct
  if (sum(abs(rowSums(A)))>.000001) print("Problem with A construction")
  if (max(diag(A))>.000001) print("Matrix A has non-negative diagonal elements")
  
  # create matrix B
  B   <- (1/dt + delta + lambda)*Diagonal(n_F) - A
  
  # initial guess for V
  V_n  <- pmax(V_low, (y + lambda * V_high) / (delta + lambda) - D_min*F_vec)
  
  # solve ODE for V recursively, starting with initial guess
  iter    <- 1
  err     <- 1
  err_vec <- c(1)
  while ((err > tol)&(iter <= max_iter)){
    V           <- V_n
    e           <- V_n-rep(V_low,n_F)
    integ_dummy <- apply(as.matrix(1:n_F),1, FUN = function(x) return(sum(pmax(V[x],x_grid)*g_vec)))
    payoff      <- y - (kappa + m)*F_vec + dx*lambda*integ_dummy
    
    # frame this as an LCP problem -- iteration n+1 as a function of iteration n
    # e(n+1) := V(n+1) - V_low >= 0 --- always the option to default and receive V_low
    # B(n) e(n+1) + q(n)       >= 0 --- HJB in continuation region
    # q(n) = - flow payoff(n) - V(n) / dt + B(n) V_low
    # e(n+1)'(B(n)e(n+1) + q(n))  = 0 and minimize quadratic form
    q       <- - payoff - V/dt + c(array(B %*% rep(V_low,n_F)))
    
    # compute e(n+1) using LCP solver
    lower_bd  <- rep(0,n_F)
    upper_bd  <- rep(Inf,n_F)
    e_n       <- LCP_fn(B,q,lower_bd,upper_bd,e,0)
    V_n       <- e_n + rep(V_low,n_F)
    
    # keep track of LCP algo error
    LCP_error     <- abs(e_n*c(array(B %*% e_n + q)))
    
    # compute error
    err     <- max(abs(V_n-V)/dt)
    err_vec <- c(err_vec, err)
    
    # update counter
    iter    <- iter + 1
  }
  
  # compute autarky value (at F = 0)
  iter2   <- 1
  err2    <- 1
  V_aut_n <- (y + lambda * V_high) / (delta + lambda)
  while ((err2 > tol)&(iter2 <= max_iter)){
    V_aut   <- V_aut_n
    G_V_aut <- G_out(V_aut, param)
    V_aut_n <- (y + dx*lambda*sum(x_grid[x_grid>=V_aut]*g_out(x_grid[x_grid>=V_aut],param)))/(delta + lambda*(1 - G_V_aut))
    err2    <- max(abs(V_aut_n-V_aut)/dt)
    iter2   <- iter2 + 1
  }
  
  # find default boundary
  def_id  <- which(V_n - V_low <= tol)
  i_def   <- min(def_id)
  F_def   <- F_vec[i_def]
  
  # compute bond price
  D <- -(V_n[(3:n_F)] - V_n[(1:(n_F-2))])/(2*dF)
  D <- c(D[1],D,D[length(D)])
  
  # compute issuance policy
  Iss <- (r-delta)*(V_n[(3:n_F)] - V_n[(1:(n_F-2))])*dF/
    pmax(conv_min,(V_n[(3:n_F)] + V_n[(1:(n_F-2))] - 2*V_n[(2:(n_F-1))]))/2
  Iss <- c(Iss[1],Iss,Iss[length(Iss)])
  Iss[i_def:n_F] <- 0
  
  # compute net issuance policy
  Net_Iss            <- Iss - m*F_vec
  Net_Iss[i_def:n_F] <- 0
  
  # create output
  output <- list()
  output[['F']]       <- F_vec
  output[['D']]       <- D
  output[['V']]       <- V_n
  output[['Iss']]     <- Iss
  output[['Net_Iss']] <- Net_Iss
  output[['V_aut']]   <- V_aut_n
  output[['F_def']]   <- F_def
  return(output)
}

# True Amador & Aguiar Model -- savings equilibrium
savings_eq_AA_fn <- function(param, num_param){
  # model such that option to default is always V_low
  # with Poisson arrival time, default value is V_high
  # model computes savings equilibrium of AA
  
  # load parameters
  V_high <- param['V_high']
  V_low  <- param['V_low']
  y      <- param['y']
  delta  <- param['delta']
  lambda <- param['lambda']
  r      <- param['r']
  m      <- param['m']
  kappa  <- param['kappa']
  
  # load numerical parameters
  n_F    <- num_param['n_F']
  F_max  <- num_param['F_max']
  F_min  <- num_param['F_min']
  
  # compute asymptotic value of V as F = F_max
  bar_D    <- (kappa + m)/(r + m + lambda)
  F_S_high <- (1/bar_D)*(y + lambda*V_high - (delta + lambda)*V_low)/(r + lambda)
  F_S_low  <- (y - delta * V_high)/r
  
  # construct grid
  F_vec    <- seq(from = 0, to = F_max, length.out = n_F)
  
  # construct value function upon jumping to F_S_high
  V_jump   <- pmax(V_low,V_low + (F_S_high - F_vec)*bar_D)
  
  # construct value function upon jumping to F_S_low or acting smoothly
  V_smooth <- V_high + (F_S_low - F_vec)
  V_smooth[F_vec >= F_S_low] <- (y + lambda * V_high)/(delta + lambda) - ((r + m)/(delta + m + lambda))*F_vec[F_vec >= F_S_low] + 
    m * (delta + lambda - r) * F_S_low / ((delta+lambda)*(delta+m+lambda))*(F_vec[F_vec >= F_S_low]/F_S_low)^(-(delta+lambda)/m)
  V_smooth <- pmax(V_smooth, V_low)
  
  # construct bond price
  D_smooth <- rep(1,n_F)
  D_smooth[F_vec >= F_S_low] <- ((r + m)/(delta + m + lambda)) + 
    (delta + lambda - r) / (delta+m+lambda) * (F_vec[F_vec >= F_S_low]/F_S_low)^(-(delta+lambda)/m-1)
  
  # construct drift rate when smooth policy
  drift_smooth <- rep(0,n_F)
  drift_smooth[F_vec >= F_S_low] <- m*(delta-r)*(r+m)*F_S_low/((delta+m+lambda)*(delta+lambda-r)) * (F_vec[F_vec >= F_S_low]/F_S_low)^((delta+m+lambda)/m+1) - 
    (r+m+lambda)/(delta+m+lambda)*m*F_vec[F_vec >= F_S_low]
  
  # solve for F_S_hat
  temp_F_I <- function(F_I){
    v_smooth <- (y + lambda * V_high)/(delta + lambda) - ((r + m)/(delta + m + lambda))*F_I + 
      m * (delta + lambda - r) * F_S_low / ((delta+lambda)*(delta+m+lambda))*(F_I/F_S_low)^(-(delta+lambda)/m)
    v_jump   <- V_low + (F_S_high - F_I)*bar_D
    return(v_smooth-v_jump)
  }
  
  # find F_S_hat
  if (min(V_smooth - V_jump)<0) F_S_hat <- uniroot(temp_F_I, lower = 0, upper = F_S_high, maxiter = 1000)$root else F_S_hat <- F_S_high
  
  # compute V across state space
  V                 <- rep(NA,n_F)
  V[F_vec<=F_S_hat] <- V_smooth[F_vec<=F_S_hat]
  V[F_vec>F_S_hat]  <- V_jump[F_vec>F_S_hat]
  
  # compute D across state space
  D                 <- rep(NA,n_F)
  D[F_vec<=F_S_hat] <- D_smooth[F_vec<=F_S_hat]
  D[F_vec>F_S_hat]  <- bar_D
  D[F_vec>=F_S_high]<- 0
  
  # compute Net_Iss across state space
  Iss               <- rep(NA,n_F)
  Net_Iss           <- rep(NA,n_F)
  Net_Iss[(F_vec<=F_S_hat)&(F_vec>=F_S_low)] <- drift_smooth[(F_vec<=F_S_hat)&(F_vec>=F_S_low)]
  Iss[(F_vec<=F_S_hat)&(F_vec>=F_S_low)]     <- Net_Iss[(F_vec<=F_S_hat)&(F_vec>=F_S_low)] + 
    m * F_vec[(F_vec<=F_S_hat)&(F_vec>=F_S_low)]
  
  # create output
  output               <- list()
  output[['F']]        <- F_vec
  output[['D']]        <- D
  output[['V']]        <- V
  output[['Iss']]      <- Iss
  output[['Net_Iss']]  <- Net_Iss
  output[['F_def']]    <- F_S_high
  output[['F_S_low']]  <- F_S_low
  output[['F_S_hat']]  <- F_S_hat
  return(output)
}

#==================================
# Value function with risk aversion
# no commitment
#==================================

risk_aversion_fn <- function(param, num_param, policies){
  # compute value function for case where government is risk averse
  # optional: specify issuance and default policies to use. If NA, just assume optimization 
  # format of those policies: a list(), with a vector of x's, a vector of g's, and a default boundary
  
  # setup for verification purposes
  dump_folder<- "C:/Users/ft.fi/Dropbox/Research/Sovereign Defaults and Commitment/Code/Dump"
  g_max      <- 5  # max issuance intensity -- stabilization of the algorithm
  
  # discretization points for numerical solution
  n          <- num_param['n_x']        # number of discretization points
  dt         <- num_param['dt']         # time interval for false transient method
  max_iter   <- num_param['max_iter']   # max number of iterations
  tol        <- num_param['tol']        # tolerance
  relax      <- num_param['relax']      # relaxation parameter
  print_flag <- num_param['print_flag'] # print flag for error checking purposes
  x_max      <- num_param['x_max']      # max debt to income for computational purposes
  
  # load parameters
  delta     <- param["delta"]
  kappa     <- param["kappa"]
  m         <- param["m"]
  mu        <- param["mu"]
  nu        <- param["nu"]
  sigma     <- param["sigma"]
  alpha     <- param["alpha"]
  theta     <- param["theta"]
  r         <- param["r"]
  gamma     <- param["gamma"]
  
  # compute the autarky welfare
  v_autarky <- 1/(1-gamma)/(delta - (1-gamma)*(mu - .5*gamma*sigma^2))
  
  # verify parameter restriction
  if (delta - (1-gamma)*(mu - .5*gamma*sigma^2) < 0) print('parameter restriction violated')
  
  # check whether policies are supplied; optim_flag = 0 if policies are supplied, otherwise optim_flag = 1
  optim_flag <- 1
  if (length(policies) > 1) optim_flag <- 0
  
  # set x_min so that we do not need to deal with singularity
  x_min     <- .001
  
  # compute x grid
  x         <- seq(from = x_min, to = x_max, length.out = n)
  dx        <- x[2]-x[1]
  dx2       <- dx^2
  
  # default barrier guess; if default policy is supplied, use that default policy
  x_bar     <- x_max
  if (optim_flag==0) {
    x_bar <- policies[['x_bar']]
    if (x_bar > x_max) print('issue with x_bar supplied: x_max < x_bar supplied')    
  }
  i_bar     <- which.min(abs(x-x_bar))
  x_bar     <- x[i_bar]
  
  # compute initial guess for debt price, value function and issuance policy; if issuance policy supplied, use it
  v_n       <- (1-(kappa+m)*x)^(1-gamma)/(1-gamma)/(delta - (1-gamma)*(mu - .5*gamma*sigma^2))
  d_n       <- seq(from = (kappa+m)/(r+m), to = .5*(kappa+m)/(r+m), length.out = n)
  v_prime_n <- rep(0,n)
  LCP_error <- rep(0,n)
  g_n       <- rep(0,n)
  # case where issuance policy is supplied
  if (optim_flag==0) {
    iss_fn        <- splinefun(policies[['x']], policies[['g']], method = "natural")
    g_n[x<=x_bar] <- iss_fn(x[x<=x_bar])
    g_n[x>x_bar]  <- 0
  }
  
  # initialization
  err     <- 1      # initialize error for value function
  err_t   <- c(1)   # keep track of errors
  iter    <- 1      # initialize iteration #
  
  # start algo
  while ((err>tol)&(iter<max_iter)){
    # if error greater than tolerance
    v        <- v_n
    d        <- d_n
    g        <- g_n
    
    # compute default option value for equity, and recovery value for debt at default
    v_interp <- splinefun(x, v, method = "natural")
    d_interp <- splinefun(x, d, method = "natural")
    v_def    <- alpha^(1-gamma) * v_interp(theta * x)
    d_def    <- alpha * theta * d_interp(theta * x)
    
    # right hand side for debt and value function
    cons     <- 1 + g*d - (kappa+m)*x
    payoff   <- cons^(1-gamma) / (1-gamma)
    RHS_d    <- kappa + m + d / dt
    RHS_v    <- payoff + v / dt
    
    # compute drift rate (invariant across iterations)
    drift_v  <- g - (mu + m - gamma*sigma^2)*x
    drift_d  <- g - (mu + m - sigma^2 - nu*sigma)*x
    
    # construct matrix A_v and A_d
    A_v      <- Matrix(0, nrow = n, ncol = n, sparse = TRUE)
    A_d      <- Matrix(0, nrow = n, ncol = n, sparse = TRUE)
    
    # diagonal elements for A
    locs      <- cbind(seq.int(n),seq.int(n))
    A_v[locs] <- -(abs(drift_v)/dx + (sigma*x)^2/dx2)
    A_d[locs] <- -(abs(drift_d)/dx + (sigma*x)^2/dx2)
    
    # lower diagonal elements -- probability to go to lower x
    locs      <- cbind((2:n),(1:(n-1)))
    A_v[locs] <- -pmin(0,drift_v[-1])/dx + .5*(sigma*x[-1])^2/dx2
    A_d[locs] <- -pmin(0,drift_d[-1])/dx + .5*(sigma*x[-1])^2/dx2
    
    # upper diagonal elements -- probability to go to higher x
    locs      <- cbind((1:(n-1)),(2:n))
    A_v[locs] <- pmax(0,drift_v[-n])/dx + .5*(sigma*x[-n])^2/dx2
    A_d[locs] <- pmax(0,drift_d[-n])/dx + .5*(sigma*x[-n])^2/dx2
    
    # impose reflection at boundaries
    A_v[1,1]  <- -sum(A_v[1,-1])
    A_v[n,n]  <- -sum(A_v[n,-n])
    A_d[1,1]  <- -sum(A_d[1,-1])
    A_d[n,n]  <- -sum(A_d[n,-n])
    
    # verify that A matrices are Markov matrices
    # if (sum(abs(rowSums(A_v)))>.000001) print("Problem with Av construction")
    # if (max(diag(A_v))>.000001) print("Matrix Av has non-negative diagonal elements")
    # if (sum(abs(rowSums(A_d)))>.000001) print("Problem with Ad construction")
    # if (max(diag(A_d))>.000001) print("Matrix Ad has non-negative diagonal elements")
    
    # construct B_v matrix
    B_v       <- (1/dt + delta - (1-gamma)*(mu - .5*gamma*sigma^2))*Diagonal(n) - A_v
    B_d       <- (1/dt + r + m)*Diagonal(n) - A_d
    
    # impose value matching at x > x_bar for d
    RHS_d[i_bar:n] <- d_def[i_bar:n]
    
    # solve for debt price, imposing value matching at default boundary
    B_d[(i_bar:n),]          <- Matrix(0,nrow = (n-i_bar+1), ncol = n)
    B_d[(i_bar:n),(i_bar:n)] <- Diagonal(n-i_bar+1)
    
    # compute d(n)
    d_n       <- c(array(solve(B_d,RHS_d)))
    
    if (optim_flag){
      # if optimization, use LCP solver to solve for optimal default boundary and optimal next period policy
      # frame as an LCP problem -- iteration n+1 as a function of iteration n
      # z(n+1) := v(n+1) - def_opt(n)   >= 0 --- always the option to default
      # B(n) z(n+1) + q(n)              >= 0 --- HJB in continuation region
      # q(n) = - flow payoff(n) - e(n) / dt + B(n) def_opt(n)
      # z(n+1)'(B(n)z(n+1) + q(n))  = 0 and minimize quadratic form
      q     <- - payoff - v/dt + c(array(B_v %*% v_def))
      
      # compute z(n+1) and thus e(n+1) using LCP solver
      z_guess   <- v - v_def
      lower_bd  <- rep(0,n)
      upper_bd  <- rep(Inf,n)
      z         <- LCP_fn(B_v,q,lower_bd,upper_bd,z_guess,0)
      v_n       <- z + v_def
      
      # keep track of LCP algo error
      LCP_error     <- abs(z*c(array(B_v %*% z + q)))
      
      # update x_high, and make sure x_bar is on the grid
      def_opt_id  <- which(v_n <= v_def + tol)
      if (length(def_opt_id) != 0){
        # update x_bar only if LCP algo did find area of state space with default exercise
        # otherwise don't update x_bar
        i_bar    <- min(def_opt_id)
        x_bar    <- x[i_bar]
      }
      
      # update debt issuance policy if we also optimize debt issuances
      v_prime_n[2:(n-1)] <- pmin(diff(v_n,2)/dx/2,-1E-10)
      v_prime_n[1]       <- min((v_n[2]-v_n[1])/dx,-1E-10)
      v_prime_n[n]       <- min((v_n[n]-v_n[n-1])/dx,-1E-10)
      g_new              <- pmin(g_max,(1/d_n)*((-v_prime_n/d_n)^(-1/gamma) + (kappa+m)*x - 1))
      
      # update issuance policy only if optimization over issuances too
      g_n                <- (1-relax)*g + relax*g_new
      g_n[i_bar:n]       <- 0
    } else {
      # if using suboptimal default/issuance policy
      # impose value matching at x > x_bar for v
      RHS_v[i_bar:n] <- v_def[i_bar:n]
      
      # solve for debt price, imposing value matching at default boundary
      B_v[(i_bar:n),]          <- Matrix(0,nrow = (n-i_bar+1), ncol = n)
      B_v[(i_bar:n),(i_bar:n)] <- Diagonal(n-i_bar+1)
      
      # compute d(n)
      v_n       <- c(array(solve(B_v,RHS_v)))
      
      # compute v'
      v_prime_n[2:(n-1)] <- pmin(diff(v_n,2)/dx/2,-1E-10)
      v_prime_n[1]       <- min((v_n[2]-v_n[1])/dx,-1E-10)
      v_prime_n[n]       <- min((v_n[n]-v_n[n-1])/dx,-1E-10)
    }
    
    # compute theoretical value of v_prime at default boundary
    v_prime_deft <- theta * alpha^(1-gamma) * v_prime_n[which.min(abs(x - theta*x_bar))]
    cons_deft    <- alpha * cons[which.min(abs(x - theta*x_bar))]
    
    if (print_flag) {
      # plot for verification purposes just in case
      setwd(dump_folder)
      
      # create list to be plotted
      plots                  <- list()
      plots[['value']]       <- v_n
      plots[['v_prime']]     <- v_prime_n
      plots[['debt']]        <- d_n
      plots[['issuance']]    <- g_n
      plots[['drift_d']]     <- drift_d
      plots[['drift_v']]     <- drift_v
      plots[['cons']]        <- cons
      plots[['LCP_err']]     <- LCP_error
      
      # plot all functions of interest
      for (plt in names(plots)){
        y_min <- min(plots[[plt]][1:i_bar])
        y_max <- max(plots[[plt]][1:i_bar])
        if (plt == 'v_prime'){
          y_min <- min(v_prime_deft,y_min)
          y_max <- max(v_prime_deft,y_max)
        }
        if ((plt == 'drift_d')|(plt == 'drift_v')){
          y_min <- -.2
          y_max <- .5
        }
        if ((plt == 'cons')|(plt == 'issuance')){
          y_max <- 2
        }
        png(paste(plt,"_iter_",iter,".png",sep=""))
        plot(x = x, y = plots[[plt]], bty = 'n', type = 'l', ylim = c(y_min,y_max), col = 'blue', xlab = 'debt-to-income', ylab = plt)
        if (plt == 'value') {
          lines(x = x, y = v_def, col = 'red')
          abline(h = v_autarky, col = 'green')
        }
        if (plt == 'debt'){
          lines(x = x, y = d_def, col = 'red')
        }
        
        if ((plt == 'drift_d')|(plt == 'drift_v'))   abline(h = 0, col = 'red')
        if (plt == 'v_prime')                        abline(h = v_prime_deft, col = 'red')
        if (plt == 'cons')                           abline(h = cons_deft, col = 'red')
        abline(v = x_bar, col = 'red', lty = 2)
        abline(v = theta*x_bar, col = 'green', lty = 2)
        dev.off()
      }
    }
    
    # compute error
    err      <- max(abs(v_n-v)+abs(d_n-d))/dt
    err_t    <- c(err_t,err)
    
    # update iteration #
    iter     <- iter + 1
  }
  
  # error message in case default option never exercised
  if (v_n[i_bar] > v_def[i_bar] + tol) print("default option never exercised -- x_max might be too low, or interest rate might be too low")
  
  # update consumption to income ratio
  cons_n     <- 1 + g_n*d_n - (kappa+m)*x
  
  # compute credit spreads
  cs_n          <- (kappa+m)/d_n - (r+m)
  cs_n[i_bar:n] <- 0
  
  # compute local consumption volatility
  cons_prime    <- diff(cons_n,2) / dx / 2
  cons_prime    <- c(cons_prime[1],cons_prime,cons_prime[n-2])
  cons_vol      <- sigma*(1-x*cons_prime/cons_n)
  cons_vol[i_bar:n] <- 0
  
  # compute stationary density
  drift  <- g_n - (mu + m - sigma^2)*x
  
  # determine point at which drift = 0
  x_zero_drift <- x[which.min(abs(drift))]
  
  # compute consumption multiplier
  cons_mult    <- (v[1]/v_autarky)^(1/(1-gamma))
  
  # construct A matrix
  A      <- Matrix(0, nrow = n, ncol = n, sparse = TRUE)
  
  # diagonal elements
  locs   <- cbind(seq.int(n),seq.int(n))
  A[locs]<- -(abs(drift)/dx + (sigma*x)^2/dx2)
  
  # lower diagonal elements
  locs   <- cbind((2:n),(1:(n-1)))
  A[locs]<- -pmin(0,drift[-1])/dx + .5*(sigma*x[-1])^2/dx2
  
  # upper diagonal elements
  locs   <- cbind((1:(n-1)),(2:n))
  A[locs]<- pmax(0,drift[-n])/dx + .5*(sigma*x[-n])^2/dx2
  
  # modify intensity matrix with treatment of lower and higher boundary
  A[1,1] <- -sum(A[1,-1])
  A[n,n] <- -sum(A[n,-n])
  
  # check that rowsum = 0 for A_temp
  if (sum(abs(rowSums(A)))>.000001) print("Problem with A construction")
  
  # compute number of consecutive haircuts given state x. 0 if x < x_high for example
  haircut           <- rep(0,n)
  haircut[x>=x_bar] <- 1 + floor((log(x_bar)-log(x[x>=x_bar]))/log(theta))
  
  # compute states post default and restructuring (with possible multiple consecutive restructurings)
  i_next           <- seq(n)
  i_next[x>=x_bar] <- apply(as.matrix(i_bar:n), 1, FUN = function(u) return(which.min(abs(x - x[u]*theta^haircut[u]))))
  
  # compute intervention matrix M
  M <- Diagonal(n)
  
  # update intervention matrix for default -- reinject defaulted firms at lower debt-to-EBITDA ratio
  locs    <- cbind((i_bar:n),(i_bar:n))
  M[locs] <- 0
  locs    <- cbind((i_bar:n),i_next[i_bar:n])
  M[locs] <- 1
  
  # store transpose of intervention matrix and FK matrix
  tA      <- t(A)
  tM      <- t(M)
  
  # create/store B matrix
  B       <- (1/dt)*Diagonal(n) - tM %*% tA
  
  # initial guess stationary distribution -- assume everything centered at bar x / 2
  f_n                 <- rep(0,n)
  f_n[round(i_bar/2)] <- 1
  
  # initiate loop
  iter2   <- 1
  err_f   <- 1
  err_f_t <- c(1)
  while ((err_f > tol) & (iter2 < max_iter)){
    # load previous period density
    f <- f_n
    
    if (print_flag) {
      # plot for verification purposes just in case
      setwd(dump_folder)
      png(paste("density_iter_",iter2,".png",sep=""))
      plot(x = x, y = f, bty = 'n', type = 'l', col = 'blue', xlab = 'debt-to-income', ylab = 'stationary density')
      abline(v = x_bar, col = 'red', lty = 2)
      abline(v = theta*x_bar, col = 'green', lty = 2)
      abline(v = x_zero_drift, col = 'orange', lty = 2)
      dev.off()
    }
    
    # update density from one loop to the next, and take into account interventions
    f   <- c(array(tM %*% f))
    
    # computing next step density using implicit scheme
    f_n <- c(array(solve(B,f/dt)))
    
    # normalize just in case
    f_n <- pmax(f_n,0)
    f_n <- f_n / sum(f_n)
    
    # normalize for distribution
    stat_dist <- f_n/dx
    
    # compute error
    err_f   <- max(abs(f_n-f)/dt)
    err_f_t <- c(err_f_t,err_f)
    
    # update counter
    iter2    <- iter2 + 1
  }
  
  # compute ergodic default rate
  lambda_d  <- .5*((sigma*x[i_bar-1])^2*stat_dist[i_bar-1]-(sigma*x[i_bar])^2*stat_dist[i_bar])/dx
  
  # compute stationary moments
  eq_object             <- list()
  eq_object[['x']]      <- x
  eq_object[['g']]      <- g_n
  eq_object[['net_g']]  <- g_n - m*x
  eq_object[['d']]      <- d_n
  eq_object[['d_def']]  <- d_def
  eq_object[['v_def']]  <- v_def
  eq_object[['v_prime']]<- v_prime_n
  eq_object[['cons']]   <- cons_n
  eq_object[['cons_vol']] <- cons_vol
  eq_object[['spread']] <- cs_n
  eq_object[['v']]      <- v_n
  eq_object[['drift']]  <- drift
  
  # compute all moments
  moments           <- data.frame(matrix(0, nrow = length(eq_object), ncol = 4))
  rownames(moments) <- names(eq_object)
  colnames(moments) <- c('mean', 'stdev', 'skewness', 'kurtosis')
  for (mom in names(eq_object)){
    moments[mom,'mean']     <- sum(f_n*eq_object[[mom]])
    moments[mom,'stdev']    <- sqrt(sum(f_n*eq_object[[mom]]^2) - moments[mom,'mean']^2)
    moments[mom,'skewness'] <- sum(f_n*((eq_object[[mom]]-moments[mom,'mean'])/moments[mom,'stdev'])^3)
    moments[mom,'kurtosis'] <- sum(f_n*((eq_object[[mom]]-moments[mom,'mean'])/moments[mom,'stdev'])^4)
  }
  
  # include in eq_object the stationary distribution
  eq_object[['dist']]  <- stat_dist
  
  # smooth functions for output purposes
  eq_fn        <- list()
  for (func in names(eq_object)){
    eq_fn[[func]] <- splinefun(x, eq_object[[func]], method = "natural")
  }
  
  # other moments to report
  other_mom                 <- c()
  other_mom['x_bar']        <- x_bar
  other_mom['default_rate'] <- lambda_d
  other_mom['d_recov']      <- d_n[i_bar]
  other_mom['cons_mult']    <- cons_mult
  other_mom['v_aut']        <- v_autarky
  other_mom['v_0']          <- v[1]
  
  # store algo performance metrics
  perf_list                      <- list()
  perf_list[['v_nb_iter']]       <- iter
  perf_list[['v_errors']]        <- err_t
  perf_list[['dens_nb_iter']]    <- iter2
  perf_list[['dens_errors']]     <- err_f_t
  
  # create output that will be returned by the program
  output        <- list(eq_fn, eq_object, moments, other_mom, num_param, perf_list)
  output_names  <- c("eq_functions", "eq_objects", "moments", "other_moments", "num_param", "algo_perf")
  names(output) <- output_names
  
  # export data in data table
  out_df        <- as.data.frame(matrix(0, ncol = length(eq_object), nrow = n))
  names(out_df) <- names(eq_object)
  for (dummy in names(out_df)){
    out_df[,dummy] <- eq_object[[dummy]]
  }
  
  # export all policies and stationary distribution into csv file
  write.table(out_df, file = paste("policies.csv",sep=""), eol = "\n", na = "NA", fileEncoding = "", row.names = FALSE)
  write.table(moments, file = paste("moments.csv",sep=""), eol = "\n", na = "NA", fileEncoding = "")
  write.table(other_mom, file = paste("other_moments.csv",sep=""), eol = "\n", na = "NA", fileEncoding = "", col.names = F)
  return(output)
}

eq_plot_fn <- function(param, eq_sol, plot_param, num_param){
  # function that plots all important functions for a given equilibrium
  # takes as given "eq_sol", output of the equilibrium function 
  # can specify plotting parameters in plot_param
  
  # download plotting parameters
  n         <- as.integer(plot_param[['n']])
  scen_name <- plot_param[['scen_name']]
  output_dir<- plot_param[['output_dir']]
  ln_thick  <- as.double(plot_param[['line_thickness']])
  ft_size   <- as.double(plot_param[['font_size']])
  use_num   <- plot_param[['use_num_sol']]
  
  # load parameters
  delta     <- param["delta"]
  kappa     <- param["kappa"]
  m         <- param["m"]
  mu        <- param["mu"]
  nu        <- param["nu"]
  sigma     <- param["sigma"]
  alpha     <- param["alpha"]
  theta     <- param["theta"]
  r         <- param["r"]
  gamma     <- param["gamma"]
  
  # download numerical solution parameter
  x_max     <- num_param['x_max']      # max debt to income for computational purposes
  
  # compute maximum debt price
  d_max     <- (kappa+m)/(r+m)
  
  # compute the autarky welfare
  v_autarky <- 1/(1-gamma)/(delta - (1-gamma)*(mu - .5*gamma*sigma^2))
  
  # download default boundary
  x_bar     <- eq_sol[['other_moments']][['x_bar']]
  
  # plot list
  plot_list          <- list()
  
  # functions to plot
  fn_to_plot <- c('d', 'v', 'spread', 'g', 'net_g', 'cons', 'drift', 'v_prime', 'cons_vol')
  plot_names <- c("debt price d(x)", "(scaled) value function v(x)", "credit spread cs(x)", "debt issuance rate g(x)", 
                  "net debt issuance rate g(x) - mx", "consumption-to-output c(x)", "drift", "v'(x)", "local (lognormal) consumption volatility")
  
  if (use_num){
    # case where we use our numerical solution calculations
    x          <- eq_sol[['eq_objects']][['x']]
    for (plt in fn_to_plot){
      plot_list[[plt]] <- eq_sol[['eq_objects']][[plt]]
    }
    # compute density
    dens_vec   <- eq_sol[['eq_objects']][['dist']]
  } else {
    # construct grid
    x_min      <- 0
    x          <- seq(from = x_min, to = x_bar, length.out = n)
    for (plt in fn_to_plot){
      plot_list[[plt]] <- eq_sol[['eq_functions']][[plt]](x)
    }
    # compute density and CDF
    dens_vec   <- eq_sol[['eq_functions']][['dist']](x)
  }
  
  # compute various outcome variables for the risk-neutral case
  risk_neutral_fn             <- list()
  risk_neutral_fn[['d']]      <- no_commit_debt_px_fn(param, x = x)
  risk_neutral_fn[['v']]      <- no_commit_value_fn(param, x = x)
  risk_neutral_fn[['g']]      <- no_commit_g_fn(param, x = x)
  risk_neutral_fn[['cons']]   <- no_commit_c_fn(param, x = x)
  risk_neutral_fn[['spread']] <- no_commit_credit_spread_fn(param, x = x)
  risk_neutral_fn[['net_g']]  <- risk_neutral_fn[['g']] - param['m']*x
  risk_neutral_fn[['drift']]  <- no_commit_drift_fn(param, x = x)
  
  # name plot names
  names(plot_names) <- fn_to_plot
  
  # name files
  file_names       <- paste(gsub(" ", "_", fn_to_plot),"_",scen_name,".pdf", sep="")
  names(file_names)<- fn_to_plot
  
  # compute y_max for density plotting purposes
  y_dens_max <- max(dens_vec)
  
  # plots
  setwd(output_dir)
  for (plt in fn_to_plot){
    # compute ylim
    y_min <- min(plot_list[[plt]][x<=x_bar])
    y_max <- max(plot_list[[plt]][x<=x_bar])
    if (plt == "cons_vol") y_max <- 2*sigma
    
    # plot all plots
    pdf(file_names[plt])
    par(mar=c(4, 4, 4, 4), mgp=c(3, 0.75, 0))
    plot(x = x[x < x_bar], y = dens_vec[x < x_bar], axes=F, xlab="", ylab="", ylim = c(0,y_dens_max), xlim = c(0,x_bar), 
         type="l", main="", lty=1, col="gray90", lwd = ln_thick)
    polygon(c(x[x < x_bar], rev(x[x < x_bar])), 
            c(dens_vec[x < x_bar], rep(0,length(dens_vec[x < x_bar]))), col = "gray90", lty = 0)
    par(new=T)
    plot(x = x[x < x_bar], y = plot_list[[plt]][x < x_bar], bty = 'n', type = 'l',axes = 'F', main = '',  xlim = c(0,x_bar),
         ylim = c(y_min, y_max), xlab = '', ylab = '', col = 'blue', lwd = ln_thick, cex.lab = ft_size, cex.axis = ft_size)
    # lines(x = x[x < x_bar], y = risk_neutral_fn[[plt]][x < x_bar], lwd = ln_thick, col = 'pink')
    grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*ln_thick, equilogs = TRUE)
    axis(2, col="black", lwd = 1, cex.axis = ft_size)
    mtext(2,text = plot_names[plt], line = 2, cex = ft_size)
    
    u <- par("usr")
    # if plotting enterprise value or q, show the no debt value
    if (plt == "v") abline(h = v_autarky, col = 'blue', lty = 3, lwd = ln_thick)
    if ((plt == "cons")|(plt == "net_g")) abline(h = 1, col = 'blue', lty = 3, lwd = ln_thick)
    
    # if plotting debt price or credit spreads, show max debt value / min credit spread
    if (plt == "d")  abline(h = d_max, col = 'blue', lty = 3, lwd = ln_thick)
    
    # include default boundary and reinjection point
    abline(v = x_bar, col = 'purple', lty = 2, lwd = ln_thick)
    abline(v = theta*x_bar, col = 'orange', lty = 3, lwd = ln_thick)
    axis(1, col="black", lwd = 1, cex.axis = ft_size)
    mtext(1,text="debt-to-income",line=2, cex = ft_size)
    box(which = "plot", bty = "l", lwd = ln_thick)
    dev.off()
  }
}

